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Synopsis 


Josephson Junction Arrays (JJAs) have been the focus of growing interest among the 
scientific community. This is in part due to a development in the standard submicron 
technology which have made the manufacture of arrays with controlled variation of 
the junction parameters feasible. Theoretically, they act as laboratory representatives 
of the much studied 2-D XY Model, which undergoes a Kosterlitz— Thouless phase 
transition. Recently the dynamics of such arrays have been studied both numerically 
and experimentally. The purpose of this thesis is to explore the dynamics of JJAs in 
detail. 

The relevant model used to numerically simulate the dynamics of JJAs is called the 
Resistively and Capacitively Shunted Josephson (RCSJ) model. In this model, the total 
current flowing through a junction is divided into capacitive, resistive, superconducting 
and noise terms. Using total current conservation at each node of a iVa, x Ny(= N) 2D 
array and neglecting magnetic shielding effects, this leads to the evolution equation of 
the phases (0) which, for overdamped junctions (zero capacitance), can be conveniently 
written in matrix notation as 

Go-' [«] = H 

Here Gq ^ is the N N connectivity matrix relevant to the system and [d] is the 
divergence vector. 

The thesis is divided into six chapters. In Chapter 1, we review the existing 
literature and place our work in its context. We develop a range of fast algorithms 
with realistic boundary conditions applicable to experimentally considered geometries 
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in Chapter 2. A general procedure consists of determining the eigenvalues and eigen- 
functions of the corresponding Green’s function satisfying the boundary conditions 
which arise out of finiteness of the array or a specific choice of the gauge (as in elec- 
trostatics). We deal with four cases and show that the complexity in the integration 
at each time step using RCSJ Model is reduced from 0{N^) to 0{N In N) for an array 
with size x Ny = N. Finally, we show that even with in the presence of defects, 
caused by bond-removal, the complexity increases marginally. 

The extension of the above procedure for the case of triangular arrays is done in 
Chapter 3. We note the equivalence of the josephson array with a lattice formed by 
sampling data at a finite number points on a 2D plane as in Digital Signal Processing. 
An extension of this idea reveals that triangular JJAs can possess two types of peri- 
odicity viz., hexagonal and rectangular. We, accordingly develop fast procedures for 
both in addition to that of the finite case. 

Chapter 4 deals with our numerical studies on defects in square arrays using the fast 
algorithms. For a linear bond disordered defect we show that vortex pinning enhances 
the system’s stability A similar behaviour is observed for the case of bond diluted linear 
defects in a long narrow array. However, the defect size is much smaller than earlier 
estimates due to finite size effects which we accordingly investigate. The energy of 
the system in the time independent regime is conveniently divided into various sectors, 
with each sector being characterised by the number of pinned vortices present in the 
array. Kink discontinuities are seen at the edges of each sector. The pinned vortices 
also cause hysteresis in the array which we examine in detail. In the time dependent 
periodic regime, waves are seen to propel these vortices forward. An investigation into 
the behaviour of distributed defects reveals that streamlined of defects have a higher 
stability as in hydrodynamics. 

We examine the case of ballistic vortex motion predicted for to occur in capacitive 
arrays in Chapter 5. Vortices, in this case, are created from a linear defect. We 
show that in contrast to a single vortex, a “vortex-street” can move through a field 
depleted region. We account for the vortex-vortex interactions including the infinite 
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array of images present across the finite array and show that the depinning current 
of a vortex is the same as in the infinite array limit. The inclusion of this term 
into the phenomenological equation of motion correctly provides the displacement of a 
single vortex as a function of time. The vortices are seen to be massless for Stewart- 
McCumber parameter < 30. Chapter 6 contains our results and discussions. 
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Chapter 1 


Introduction 


Superconductivity has assumed an important role in today’s technical world with its 
usage in the manufacture of scientific devices like superconducting magnets and wires, 
SQUIDS, etc. In this regard, granular superconductors, in which micrometer sized 
superconducting grains are placed in a non-superconducting host, hold special impor- 
tance. Apart from their natural occurrence in the manufacture of bulk superconductors, 
they are now routinely fabricated in various laboratories using highly sophisticated 
submicron technology. Their two-dimensional counterparts consist of a network of 
Josephson links connecting nearest neighbour superconducting grains[l]. Its industrial 
applications apart, the study of such systems is important from both theoretical and 
experimental viewpoints. Recent interest in the study of such arrays owes much to the 
development of standard submicron technology which permits one to fabricate such 
arrays with relative ease as well as vary their parameters over a wide range of values in 
a controlled manner. Moreover, these systems are experimental realisations of theoret- 
ically well-studied models and provide much insight into the behaviour of conventional 
and granular high-Tc superconductors. 
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1.1 Properties 

Since a Josephson link consists of superconductors grown on a host material, they are 
classified into three categories in accordance with the material used as the host. These 
are : 

SIS: Here an insulating material (e.g. Si ) is used as a host while the superconductors 
are made of Nb. 

SNS: The superconducting islands (e.g. Pb) are placed on a host which is a normal 
metal (e.g. Cu). 

WL: A weak link consists of a narrow superconducting bridge between two supercon- 
ducting nodes. 

The above classification holds for Josephson Junction Arrays (hereafter referred to 
as JJAs) too, since they are created from josephson links. All the categories of arrays 
listed above are operated well below the BCS temperature of the superconducting 
islands. The properties of the SNS and the SIS type JJAs are very similar and these 
are what we would be mostly discussing here. 

The procedure used for fabricating such arrays[2] consists of creating a stencil or 
mask using photo or electron beam lithography and the subsequent etching of the 
photoresist material. The required material, e. g. the superconductor or the insulator, 
is next deposited through this mask by evaporation followed by a “lift off” of the stencil 
layer. Various experimental facilities have used different combinations of materials 
to achieve the desired properties of arrays and we list some of the basic properties 
in Table. 1.1. Usually square arrays are manufactured. Recently, however, the Delft 
Technical University, Netherlands [3] and others[4] have manufactured triangular arrays. 
Other exotic types of geometries e. g. Penrose tiling schemes, Amman lattices and 
Strepenski gasket lattices have also been designed and studied[5]. To describe the 
system at hand, we label each superconducting site (also called a node) by f, a two 
dimensional vector (see Fig. 1.1). Noting that f can take on discrete values as specified 
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Figure 1.1: A square lattice with the light shaded regions depicting the superconductors 
and the dark shaded ones the insulator, r indicates the sites while R represents the 
dual sites. 

by the lattice, we also interchangibly use the scalar index i to describe the node, since 
there exists a definite relationship between the two indices, specified by the underlying 
lattice structure. 

The essential macroscopic physics of such arrays is captured if one considers that, at 
temperatures below the BCS temperature, the wavefunction of the condensed Cooper- 
pairs at the superconducting sites, ^(r, t), to be of the form 

= \'^\expi9(r,t) (1-1) 

This is a oversimplification as it ignores the variation of |'^'| with t. Moreover, it 
inherently implies that the phase 0{f, t) is the same over the entire superconducting 
island. However, as shown by Eikmans[6] the phase variation over a typical island 
(denoted by Rq in Table. 1.1) is ~ 10“'*. Under these approximations, the whole array 
is described by a single parameter set {^(F, t)}. 

Having laid down the basic tools for describing the system under consideration, we 
note that the study of such arrays can be broadly classified into two streams. The first 
stream consists of studying the thermodynamic properties of arrays in statistical equi- 
librium i.e. the static properties, while the other aims at investigating the behaviour 
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of JJAs far from equilibrium i.e. the dynamical properties. We discuss both directions 
briefly. A comprehensive account, including a microscopic treatment of such arrays is 
present in books of superconductivity [7, 8] and standard text books on JJAs[9]. 

The Hamiltonian governing the statics of such arrays is; 


= Jo cos(6{f, t) - 6{f\ t)) (1.2) 

{ff'> 

where {f, r') represents the nearest neighbours of f and Jo = (^*o)/(2e). Here Zq is the 
single junction critical current. 

One assumes that iq « i^o where Zco is the critical current for which superconduc- 
tivity is destroyed within the islands. This holds in the weak coupling limit for arrays, 
which in turn means that the thickness of the barrier between the superconducting 
sites is ~ 1°A. To keep the nodes from becoming normal, the arrays are operated well 
below the BCS temperature. 

Perhaps the most fundamental relation, in the context of JJAs, is that between the 
voltage (AV) developed across a junction due to the existence of a phase difference 
(A^) at the the two superconducting sites of a junction. This is called the Josephson 
relation [1] and is given by 


Ay = 


h dM 
2e dt 


(1.3) 


The relationship has found an important application in the setting up of voltage stan- 
dards since it involves fundamental constants. 


1.2 Arrays In Thermal Equilibrium 

Knowing the Hamiltonian (Eq.[1.2]), the thermodynamic properties of the system can 
be calculated from the canonical partition function Z 

S = £cWexp(-^) (1.4) 

where the D{6) represents integration over all the field variables {0(i^} and Ub is the 
Boltzmann factor. 
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Using generalised tlieorems[10-12], forbidding any long range order (LRO) for di- 
mensions D < Dc = 2, it was predicted that JJAs (as well as helium and super-fluid 
films) cannot undergo any phase transitions. The theorems suggest that for D < Dc, 
the fluctuations in the system destroy LRO. The absence of LRO rules out the possibil- 
ity of an ordered phase. This is, however, not the case. Berezenskii[13] and Kosterlitz 
and Thouless[14, 15](BKT) showed that in such continuous-sjunmetry systems a tran- 
sition of the “topological excitations” or vortices can occur: At low temperature the 
vortices and anti-vortices present in the system are bound together to form dipole pairs 
which dissociate at high temperatures. To understand the new physics involved, we 
digress for a moment to discuss the new concept of vortices. 

Physically, a vortex signifies a region of circulation of the concerned field variables. 
In hydrodynamics the vorticity u is related to the velocity fields u as l 3 = V x u. 
However, in the case of JJAs, there exists a quantisation condition describing a vortex. 
Mathematically, as a result, vortices are defined by a constraint in the phase differences 
of the field variables, 9i — 6j. Any path enclosing the center of a vortex is constrained 
by 

^ Oi — 6j = ±27r7n where m is an integer (1-5) 

where the phase differences 6i — 9j have range — tt < 9i — 9j <7r. Since the constraint 
is over any path, the vortices are seen to be extended objects. A similar condition is 
used to describe a vortex in helium films where the phase angles 5(f) = expi9{r) are 
quantised as ^'V 9 ■ dl = 2rmr. 

By mapping the hamiltonian for JJAs onto the 2-D XY model [16] one can extract 
the critical properties of the arrays. Many authors have performed detailed calculations 
in this regard[17-21] and we present a brief overview of the procedures followed. 

Using the dual-transformation approach[19, 20], one can factor Z as 

2 = ( 1 . 6 ) 

where Z^^ is the spin— wave part and Z^^ the Coloumb gas part of the partition 
function. The spin— wave part is irrelevant to the study of critical behaviour while 
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Z^^=y:exp irKoTm{7^mif')ln (1-7) 

L J 

where «A"o = tJoHUbT), oq is the lattice constant of the array and yo, the bare fugacity 
of the system, is given by 

Vo = exp(-jr=jr„/2) (1.8) 

In arriving at Eq.(1.7) one uses the neutrality condition 

2m(f) = 0 (1.9) 

r 

which arises from the tf(l) symmetry 6i —^di + a{t) of H (Eq.[1.2]). 

Using a renormalisation technique[14] one can then derive the renormalised cou- 
plings (Ki) and fugacities (y^) at an arbitrary length scale I in the form of differential 
equations (there is a third differential equation for the free energy of the system). The 
linearisation of these differential equations provides the renormalisation flows and the 
transition temperature (Tbkt)- H can be shown using this procedure that the scaled 
dielectric constant K^o/Kq attains a universal value of 2/7r at Tbkt before dropping to 
zero and that the phase correlation has a power law behaviour, ~ r’’ where 

v(Tbkt) = 1/4. 

Experimentally, Sanchez etal.[22] observed a square root dependence of temper- 
ature of the dc sheet resistance indicative of a BKT transition. Since then many 
workers[23-30] have seen signatures of the BKT transition in such arrays in the pres- 
ence/ absence of an external magnetic field. Another method to observe the transition 
is based on the fact that vortex— antivortex pairs experience a Lorentz force in opposite 
directions in the presence of an external current drive. This causes the pairs to disso- 
ciate and flow at right angles to the current. This in turn results in a voltage drop V 
across the array which is found to be related to the input current / by 

V ~ a(T) = 1 4- 7rK^(T) (1.10) 

Since K^o has a universal drop at Tbkt, a(Tgj^-j-) = 3 and a(r^;^.j,) = 1[31]. Both 
experimental studies[32] and Monte— Carlo (MC) simulations[30] based on dynamic 
theories[33] have confirmed the above universal behaviour. 



CHAPTER 1. INTRODUCTION 


7 


1.3 Imposition of Magnetic Field 

When an external magnetic field is applied to a JJA, the line integral A(r, r') of the 
vector potential A enters the hamiltonian through the minimal coupling prescription 
{0(f) — 9{r')} —* {0(f) — 0(f*O ~ 

H — cos{0(f',t) - 0(f",t) - A(f, f')} (1-11) 

(r"') 

where j4(f', r') = Sr ^ ■ dl. Here, (f)Q = (hc)j(2e) is the elementary flux quantum. 
The magnetic field penetration depth of the superconducting islands is assumed to be 
much smaller than the size of the islands, which are thus viewed as completely expelling 
the magnetic field. Moreover, the flux piercing the experimentally prepared junctions 
is taken to be much smaller than 0o- This ensures that the magnetic shielding effects, 
arising from the variation of the phase around the superconductor-insulator contact 
plane, become negligible. These effects lead to a suppression in io- The vorticity 
condition in the presence of magnetic fields is 

E {6i-dj-Aij) = 2Trf (1.12) 

where / is the flux per plaquette in units of ^o- In both Monte Carlo (MC) and 
dynamical simulations of JJAs in the presence of a perpendicular magnetic field, the 
vector potential is incorporated into the governing equations using the Landau gauge 
condition 

{ 27r/ for horizontal bonds 

(1.13) 

0 for vertical bonds 

Following an approach similar to the one outlined above (Sec. 1.2) one can write 
down 2 as in Eq.[1.7] by replacing m(f) — fl(f) + / where Q.(r} are integers and / is 
the static contribution by the magnetic field. The neutrality condition also reads 


Sn(f) + / = o 


( 1 . 14 ) 




Moreover, since Z has the symmetries 


f f + n n E I and / — + — / 

the only relevant range of / is [0, 1/2]. In fact, the symmetry of / over integral values 
is often used to test the quality of the experimentally prepared films[34]. 

The fact that magnetic fields induce frustration[35, 36] in an array is most easily 
seen for the case of / = 1/2. Here alternate horizontal bonds have A{j — it (and 
7i = —Jocosdi — 9j) while the rest have Aij = 0 (and 7{ = Jq cos Oi — 9j). Thus each 
plaquette has three anti-ferromagnetic and one ferromagnetic bond causing circum- 
stance frustration. The case of / = 1/2 is called the fully frustrated case. The ground 
state configuration for such an array, even for T = 0, is much more complex since 
field-induced vortices are created to satisfy the neutrality condition (Eq.(1.14)). Teitel 
and Jayaprakash[37], in their pioneering work, showed that iox f = p/q [p and q are 
integers and coprime) the ground state isaq-xq periodic structure of the vortices. The 
upper bound of the critical current passing through the junction was estimated to be 



This implies that for T = 0, ic — s- 0 independent of the size of the system. Halsey[38] 
considered the “staircase” state of such d^qx q lattice and showed that ic{f = pi q) = 
ic{f = 0)q~^ independent ofp. The “staircase” state is, however, not the general ground 
state as revealed by Monte Carlo studies conducted by Teitel[39]. For irrational p/q, 
it has been predicted that the ground state configuration is that of a superconducting 
glass [40H12]. 

The particular case of full frustration where / = 1/2 has received a great deal 
of attention as a realisation of Villian’s “odd-model” [35]. Although the relevant 
theory[43-47] and experiments[16, 26, 29, 48] are still far from complete in this area. 
The ground-state configuration, in this case, consists of “chequerboard” pattern of 
±1/2 vortices. The existence of an unnatural phase transition - a mix of the BKT and 
the Ising phases, as also the universality class of such a transition have been points of 
considerable debate. 
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1.4 Capacitive Arrays 


Arrays (especially SIS arrays) can have appreciable junction capacitance Cj as well as 
junction to ground capacitance Cq- The corresponding charging energies for Cooper 
pair transfer are Ecg = (2e^)/C'(3 S'lid Ecj = {2e^)/Cj respectively. Such capaci- 
tances tend to destroy phase coherence in the ordered phase. In fact for Ecg ^ Ej, 
the junction coupling strength, quantum fluctuations are expected to predominate. 
Theoretically, the treatment of capacitive arrays follows the route of writing down an 
effective action[49] 


S_ 

h Jo h 



( 1 . 16 ) 


where = l/(kBT), m = hJ /{2e^)CG^ M = ftJ /{2e^)Cj and r is discretised in units of 
tq. One then has an N x N x Nr lattice with Nr = ^hfro. Thus capacitive arrays in 
d dimensions are effectively treated by a d -I- 1 dimensional field theory. As noted by 
Anderson and Abeles[50] such arrays can have a Mott-like insulating phase at T = 0. 
Inclusion of perturbative effects is seen to bring down Tba't[ 51]. In fact, such arrays 
are seen to undergo a novel first order Quantum Induced Transition for temperatures 
well below Tbkt[^‘J\- 

The dynamics of vortices in such arrays has been the focus of much recent interest 
It has been speculated by Eckern[53, 54] and Korshunov[55] and that vortices in such 
arrays have a finite mass and under certain conditions can move ballistically in a field 
depleted region. Van der Zant et al.[3] claim to have seen ballistic motion in triangular 
capacitive arrays. Numerical simulation [5 6], however, does not show any such motion. 
Although various reasons have been provided to account for the discrepancy by Geigen- 
muller et oZ. [57], Yu et aZ. [58, 59] and Fazio et aZ.[60] the reason for the discrepancy is 
still unclear. 
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Figure 1.2; The figure shows how the phases in the left plaquette have to rotate by a 
finite amount so that the vortex (represented by +) can move to the right plaquette. 


1.5 Dynamics of Arrays 


As we have seen previously, an application of external currents introduces the concept 
of phase slippage and a voltage drop in the array (Eq.(l.lO)). Thus the motion of a 
vortex present in such an array, causes dissipation. However, in contrast to continuum 
superconductors, where an infinitesimal external current causes motion of the vortex, 
JJAs provide a pinning potential which restricts such motion. This is due to the 
inherent discrete nature of the arrays. The vortex present at the center of a plaquette 
is in a low-energy site, and has to overcome a finite potential barrier to move to 
another plaquette (another low-energy site). To see this more clearly, we note that 
the motion of a vortex from a plaquette to the next must cause a finite rotation of 
the phases which translates to a definite energy barrier (see Fig. 1.2) opposing such 
motion. Using a pure sinusoidal potential for the pinning, Lobb et al.[61] have shown 
that vortices in an infinite array can move only if the external current scaled by the 
single junction critical current is ~ 0.1. 

To study the dynamics of JJAs, we take recourse to the well known Resistively 
and Capacitively Shunted Junction (RCSJ) Model[62]. More refined models exist[9] 
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to solve the single junction dynamics but the RCSJ model is, for most purposes, com- 
pletely adequate. According to this model, the total current i(r, r ', t) between the 
superconducting island f and f ' in general can be divided into four “channels” as 


i(f, f ', t) = zs(f, f', t) -i- r\ t) 4- ic(r, r\ t) + ipir, f\ t) (1-17) 

All the currents are scaled in units of The superconducting component, is{r^r'^t) 
is related to the phases 0(r, t) by the relation 


t', t) = sin[0(f) — d{r") — A(f, r ')] 


(1.18) 


and is due to the tunneling of Cooper pairs across the junction. In the RCSJ model, we 
consider a linear relationship between the voltage developed due to the current flowing 
through the junction i.e. we shunt the junction by a ohmic normal-state resistance Rn- 
This linear approximation adequately describes the dynamics even though the shunt 
resistance at the operating temperatures (< Tbcs) is not equal to Rn. In Eq.(1.17), ic 
represents the current due to a flnite inter-junction capacitance, C. Thus 


leioRn 


{6{r) — O) 




(1.19) 

( 1 . 20 ) 


Finally, there is the fluctuation current which we take to be due to a white 

noise source, with no spatial correlations, i. e. between the noise in different junctions 


{ipir, r t)iF(Fi, - t') - f{) - 6(r - fl)6{f' - ri)] 


/3u, 


( 1 . 21 ) 


Here — (2eioRn)/S- and /? inversely proportional to the temperature. Since the 
temperature enters through (3 = J{T)/{kBT) = hio(T)/(2ekBT), Eq.(1.21) is used in 
dynamical simulations of arrays at finite temperatures. 

The time evolution of the phase difference = 9{r) — 6(r') across a single junction 
is thus governed by the equation: 
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+ A^(r) + sin(A^(f) - 27rA(f, f ')) + 2>(f, f', t) = z(r, r ') (1.22) 


Here we measure time in units of 1/wc- The Stewart-McCumber characteristic 
parameter[63, 64], ,8c is defined in terms of the frequency, Wc, of a single junction: 



and the plasma frequency 



(1.23) 


(1.24) 


In practice the characteristic frequency is of the order of One easily finds 

that Pc is the typical charging energy of the junction, C{IcRnf , divided by the typical 
coupling energy, J . For the SIS arrays. Pc can be varied continuously between /?c -C 1 
(overdamped junctions) and Pc 1 (underdamped junctions). The SNS arrays are 
always strongly overdamped. 

In arrays, these single-junction equations are coupled through the condition of 
current conservation at each superconducting node of the array. The. resulting Total 
Current Conserving (TCC) dynamics[33] leads to a set of iV, first-order, coupled, non- 
linear differential equations for the case of an array with Nx x Ny = N overdamped 
josephson junctions. We shall in Chap. 2 show how such a set of equations can be 
written such that they are amenable to theoretical analysis and simulations. 

Both experiments and numerical simulations have in tandem probed various pa- 
rameters of the array in detail. Both overdamped (/?<- = 0) and underdamped (finite 
Pc) current driven arrays have been studied at zero and finite temperatures in the 
presence/ absence of magnetic fields. We briefly discuss some directions in which inves- 
tigations have proceeded below. 


1.5.1 Disordered Arrays 

Disorder is introduced in JJAs[65] in the form of 
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• Bond disorder, where the Josephson coupling J — with J,j following a speci- 
fied distribution (usually Gaussian), 

• Bond dilution, wherein bonds are removed either randomly or systematically to 
form a defect and 

• Positional disorder, wherein the magnetic vector potential is a random or 
Gaussian variable. 

Bond disorder has been used to distinguish between a BKT and Ising like transi- 
tions[66, 67] for / = 1/2, to observe a lowering in the resistivity due to pinning[68, 69], 
to facilitate the movement of vortices[3] and to simulate spin glass like phases in 3D 
arrays[70]. Positionally disordered arrays, have been studied[65, 71, 72] to look for a 
predicted[73] reentrant glass like phase. A novel vortex pattern[74] and a plastic-like 
flow of vortices are seen in such arrays. Disorder in randomly diluted arrays have been 
seen to preserve the scale invariance of the BKT transition[75]. The breakdown of a 
current driven array in the presence of a bond diluted linear defect is seen to occur 
when a vortex path spans the entire sample[76, 77]. The system becomes chaotic as 
the “vortex-street” spreads to various columns. Chaos also occurs due to a mixing of 
positive and negative vortices[78, 79] in an array driven by a linear current profile and 
in capacitive arrays in the presence of electro-magnetic radiation [80]. 


1.5.2 Shapiro Steps 


The I-V characteristics of a single junction driven by external sinusoidal current I = 
Jr/sin27rz/t contains the the well known steps first observed by Shapiro[81]. These 
result from the fact that the voltage across the junction is mode locked at well known 
values given by Ch. Leeman et aZ.[48] and T. D. Clark[82] experimentally 

observed that for b, N x N array these steps were amplified by a factor of N and that 
the average voltage of the step is 


w„ = 


nNhv 


n = 0, 1, 2 


N = 1000 


(1,25) 


2e 
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Under the application of an external magnetic field characterised by / = p/g (p and q 
are integers and coprime), fractional Shapiro steps were observed by Benz ei al.[83] for 


(^)n = 


nNhv 

2qe 


(1.26) 


These authors explained the steps in terms of a motion of the vortex lattice as a 
whole, by an integral multiple of the lattice constant a. Numerical simulations[84] for 
f = 1/q revealed additional steps at (V)^ = theoretically shown that 

the extra shapiro steps would occur if certain rows of the vortex lattice, as opposed to 
the lattice as a whole, moved by qa. 


1,6 Overview Of Current Research 

Josephson Junction Arrays are “toy universes” which exhibit a broad range of inter- 
esting phenomena[85-87]. As a result, the study of such arrays has developed along 
various lines. We provide a brief overview of the related topics below. 

Since JJAs provide laboratory models to verify and test the predictions of the 
2D XY Model much effort has been devoted to observing the critical behaviour of 
such arrays. These experiments have thus helped to verify and extend the theoretical 
underpinnings of the field. A considerable part of the research has been devoted to 
understanding the case of full frustration [3 7, 44, 45, 47, 88-97]. Disordered samples 
have also been studied in the presence of external magnetic fields since they are believed 
to be similar to spin-glasses [40, 42, 41, 46]. Generalisations of this model have also 
been attempted[96, 98-100]. Theoretical predictions[73] of reentrant behaviour have 
yet to be substantiated unambiguously[65,.71, 72]. 

The controlled fabrication of arrays have furthered investigations of their dynam- 
ical properties[30, 68, 76, 78, 101-106]. These have resulted among other things, in 
the observance of flux-flow behaviour, vortex pinning[77, 104, 106] and more compli- 
cated vortex dynamics. All of these have direct implications on the superconducting 
properties of granular superconductors. As is the case of statics, the dynamics of these 
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“collective vortex-variables” have been correlated with breakdown phenomenon and 
with the advent of chaotic behaviour in JJAs[79, 80, 107]. It has been shown that 
for disordered samples their exists a novel pattern in which the vortices arrange them- 
selves at 27° angle to the current direction far from equilibrium [74, 108]. Moreover, 
a plastic-like flow of vortices is seen to occur in such highly disordered samples[109, 
110 ]. 

The discovery of the Josephson eff^ect has triggered a large effort at using such 
junctions as generators for microwave signals[9] with the tunable frequency related to 
the input dc drive current. However, the low power supplied by such junctions makes it 
inappropriate for industrial applications. Coupled arrays of such junctions capable of 
producing a coherent radiation was an alternative to overcome such a limitation. These 
required the existence of in-phase solutions primarily studied in 1— D arraysjlll, 112]. 
So far, no coherence has been observed in 2-D dc-driven arrays[113, 114] in contrast 
to rf-driven arrays[115, 83]. 

A regime in which the charging energies of the junctions (especially in SIS arrays) 
is comparable in magnitude to the Josephson coupling has opened up a wealth of new 
and interesting concepts and observations[116, 117, 114, 118]. Recent experiments have 
accordingly probed superconductor-insulator transition[118] previously predicted and 
observed in granular superconductors. Arahranov Casher[119] and Bohm[120] effects 
have also been seen in JJAs. 


1.7 Outline Of The Thesis 

Here we present an outline of this thesis. 

We shall be mostly concerned with the dynamics of arrays in the presence of cur- 
rent drives and defects in this thesis. Our approach will be to combine the theoretical 
analysis and numerical investigations to develop the machinery of integrating the cou- 
pled RSJ equations in Chap. 2. The machinery includes a range of fast integration 
algorithms which take into account realistic boundary conditions[121]. In Sec. 2.2 we 



present a simple derivation of the Eikmans— Himbergen Algorithm[122] using the tech- 
nique of Green’s function and cast the RSJ equations into the familiar language of 
electrostatics. We also clarify the roles played by gauge invariance and boundaries in 
the context of JJAs. These techniques are generalised to deal with the cases of single 
and double busbars in Sec. 2.3. Finally, in Sec. 2.4 we extend the scope of the fast 
procedures to that of a single broken bond and more generally to arrays with linear 
defects. 

We extend the ideas developed in Chap. 2 to design fast algorithms for triangular 
arrays[123] in Chap. 3. In Sec. 3.2 we show how the Dirac bra-ket algebra can be used to 
describe the JJA. In developing the algorithm we explicitly demonstrate the relevance 
of ideas used in digital signal processing to the systems at hand (Sec. 3.3). We show 
that such arrays admit two types of periodicity and develop fast computational schemes 
for both. Finally, in Sec. 3.4, we extend the scope of the computational procedures to 
a finite triangular array. 

In Chap. 4, we study the effect of bond disordered/ diluted JJAs[124, 125]. The 
case of a linear defect is dealt with in Sec. 4.2. It is seen that for the case of a bond 
disordered linear defect, vortex pinning enhances the stability of the system against 
breakdown from superconducting to resistive flow. In the bond diluted linear defect 
case, we observe pinning by the perfect lattice for defect sizes much smaller than earlier 
predicted. A study of the finite-size effects shows scaling behaviour. The energy of 
the array in the presence of pinned vortices is seen to undergo discontinuous jumps 
and to display hysteretic behaviour, which we investigate in detail in Sec. 4.2.3. The 
periodic regime of such arrays is shown to consist of waves which propel the vortices 
forward. These are traced in Sec. 4.2.4. Sec. 4.3 deals with the idea of streamlining a 
defect to enhance its stability of the flow as for hydrodynamics. We show that even 
with increased bond dilution one can defer the onset of the resistive transition. 

Motivated by the discrepancy reported in the literature between numerical observa- 
tions, in Chap. 5, we take up the study of ballistic vortex motion in JJAs[126]. Sec. 5.1 
presents a configuration in which vortices are created by defects rather than magnetic 
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fields. We show in Sec. 5.2 that although a single vortex does not move in a field de- 
pleted region, a “vortex-street” does so even in overdamped arrays. This is primarily 
because the vortices, seen as topological charges, interact via a long-range potential 
with one another (Sec. 5.3). We take into account the efi'ect of image charges present 
across the finite array and numerically show that the motion of a single vortex obeys 
a phenomenological equation of motion containing the vortex-vortex interaction term 
(Sec. 5.4). 


Appendix l.A 


Description 

parameter 

units 

Nb-Si-Nb 

A.l-Al0a;^“"A.l 


Tco 

K 

9.50 

1.14 


Tkt 

K 

5 

- 

Critical Current 

io 

fxA 

1 

0.5-0.005 

Normal-state Resistance 

Rn 

kQ 

0.5 


Critical Magnetic field 

H^{T = 0) 

Gauss 

1980 


Magnetic Field 


Gauss 

0.4 

0.4 

Linear size 

Is 


5 

5 


is 

— 

0.2 

0.2 

Area of single junction 

Aj 

fj,m? 

0.01-0.1 

0.01-0.1 

London pentration depth 

Xq 

fim 

0.039 

0.016 

Coherence length 

Co 


0.04 

1.60 

Variation of phase in island 

Ro 

- 

10-4 



Table 1.1: Overview of the important array parameters. 
























Chapter 2 


Fast Algorithms for Square Arrays 


2.1 Introduction 


With the size of experimental arrays increasing continuously, and with the number 
of interesting effects best seen only in large arrays going up in equal measure, it has 
become imperative to find ever more efficient algorithms for implementing the corre- 
sponding simulations, inclusive of all the experimental conditions. An example of the 
latter for current-driven arrays is the presence of bus-bars, through which the external 
current can be conveniently injected or withdrawn. 

These algorithms must address the problem of integrating the complete set of cou- 
pled equations for each superconducting site describing the array in terms of the phase 
variables, Within the framework of the RCSJ Model, the equations found in Sec. 1.5, 
have to satisfy current conservation at each site[33]. This yields the form 


E 

(•i) 


chcpeij 

2e dt^ 


+ 


h ddij 
2eR dt 


+ io sin 6ij 


= Vi 


( 2 . 1 ) 


where C, R and io are the shunt-capacitance, shunt-resistance and critical current of 
the junction respectively. We note that Eq.(2.1) holds under assumptions of zero tem- 
perature, zero magnetic field and infinite perpendicular magnetic penetration depth. 
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The last of these allows us to neglect self-induced magnetic field efliects. We can rewrite 
Eq.(2.1) as 

(ij) 

where {iext)i — is the normalised current being fed to or extracted from the array 

site 2 , 3c = 2eioR^C/h is the McCumber-Stewart parameter and t = t(2eioR)/h is the 
time measured in units of the characteristic period = hl{2eiQR). The summation 
on j, over the nearest neighbours of i, can alternatively be expressed in terms of a 
multiplication by the discrete laplacian, Gq^. Eq.(2.2) then assumes the matrix form 

= ( 2 - 3 ) 

If we set di = X, and 6i = yi, the complete set of dynamical equations reads 
.i/j = —di([x], [y]) and x,- = y,-. For the overdamped case, corresponding to 
3c = 0, the relevant equations are 

(GJ‘) = -d, (2.4) 

3 

where — d,- is now given by — S)(,j) sin but Gq^ is, of course, the same as in 
Eq.(2.3). 

It follows that for an. Ny = N array each integration time step of Eq.(2.4) has 
a complexity 0{N^). This is because at every upgradation of the N state variables, 
6 in the underdamped, and 6, in the over-damped case, the constant N x N matrix 
Go has to be multiplied by the divergence vector [d\. It was first noticed by Eikmans 
and Himbergen et al.[122] that the form of the is such that this multiplication can 
actually be carried out in 0{N In N) steps. The procedure was subsequently improved 
upon by Dominguez et al.[74] who combined the fast-fourier transformation used by 
Eikmans et a/.with the method of cyclic reduction to achieve a roughly three-fold 
increase in speed. This algorithm was also used to study JJAs with inductance [110]. 

It is noteworthy that these algorithms are applicable even in the presence of an 
external magnetic field. Indeed, the application of a field, Bqz, perpendicular to the 
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array transforms the phases dij into the gauge-invariant combinations 6ij — 2TrAij where 
A,j = l/<pol2{ A - dl, A is the vector potential and (po = ?t/(2e) is the elementary flux 
quantum threading a plaquette. This transformation clearly affects only the divergence 
term, and leaves untouched the matrix , on whose form the algorithms are based. 
Similarly, white noise, which is taken into account by introducing a noise current into 
Eq.(2.3), also modifies only the divergence and hence does not affect the applicability of 
these algorithms. The statement continues to hold even if we introduce bond disorder, 
i.e. make io junction-dependent (see Chap. 4). 

In this chapter we extend these fast algorithms firstly to arrays with busbars and 
secondly to those with defects in the form of missing bonds. Busbars often form 
the current injection and/or extraction edges of experimental arrays. In experiments 
involving vortices, which are repelled by busbars, these have been used to produce 
collimated vortex-streets[3]. The dynamics of vortices have also been investigated with 
one edge shorted by a single busbar[78, 79]. Arrays with defects have likewise arisen 
in a number of contexts. The breakdown of superconductive flow in current-driven 
arrays with linear defects, for example, has been investigated in some detail. The 
exploration of the multiple-vortex sector, which arises in this study, requires running 
on large arrays and is all but impossible without the algorithm we discuss in Chap. 4. 
Defects can also be used to provide a collimated beam of vortices (see Chap. 5). 

The chapter is organised as follows. In Sec. 2.2, we present a simplified derivation 
of the Eikmans-Himbergen algorithm and interpret all the key equations in the famil- 
iar language of electrodynamics. Apart from being more transparent, our derivation 
clarifies the role of boundaries and the connection between the lattice and continuum 
descriptions of the systems being studied. These insights are used in Sec. 2.3 to gen- 
eralize this algorithm to arrays with busbars. The case of single busbar forces us to 
resort to the technique of cyclic reduction, which we consequently discuss at this point. 
In Sec. 2.4, we extend these algorithms to JJAs with defects, created by eliminating or 
adding bonds. Sec. 2.5 contains a summary and discussion of our results. 
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2-2 The Eikmans— Himbergen algorithm 

2.2.1 Preliminaries 

All efficient algorithms for the numerical integration of Eq.(2.4) make essential use 
of the properties of the discrete laplacian, We thus begin by listing the more 

important properties of this matrix. From its definition (as given by Eqs. (2.3 and 2.4), 
it follows that Gq ^ specifies the connectivity between different sites. More precisely, 
= (^ 0 ^)^-. = i (site i connected to site j) and 0 otherwise. The 

diagonal element ^Gq^). = total number of points to which the site i is connected, 
and hence rr(Go^) is twice the number of bonds in the network. 

Clearly, Gq ^ is real and symmetric. As a result, its eigenvalues are real and its 
eigenbasis complete in the iV-dimensional space of states. Furthermore, det(Go = 0. 
This can be deduced as follows. We first note that since Eq.(2.2) involves only the 
phase differences (and their time derivatives) it is invariant with respect to the global 
transformation 9i — + q:(t) V i. Substituting this into Eq.(2.4), we have 

E + “(’■)) = -i. = E (Go-),/,- 

i j 

This implies that d:(r) (Gq^].. = 0 V i. I.e., if we fix i and sum over all the 

columns we get zero. (The same is true, of course, if we fix j and sum over i since 
the matrix is symmetric). Hence det(Go^) = 0 and only {N — 1) of the N equations 
in Eq.(2.2) are independent. (We could, alternatively, note that ('Fo)i = 1 V i is an 
eigenvector of Gq ^ with eigenvalue 0). 

To explicitly evaluate 0,- from 6i we have to eliminate the extra degree of freedom. 
For what immediately follows, the most convenient choice is J2i — 0- Eq.(2.4) is then 
replaced by 

E(S„-'),/i = -Bi (2.5) 

3 

where D,- = d,-, i = 1 • • • (TV — 1) and Dn = 0 while i = l”'(W — 

1), j = 1 N and . = 1 V y. Moreover, the vector [d] must satisfy d,- = 0 as 
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can be seen by performing an additional sum over i in Eq.(2.4) and using the fact that 
E, (Go ^),-j = 0. For JJAs this condition is automatically satisfied since % and smd,j 
are odd functions of their arguments and the net external current fed to the array is 
zero. 

One can, in principle, invert and determine 6 from [Z)]. This is the usual 
0{N‘^) process. We shall refrain from adopting this procedure and return, instead, to 
working with Eq.(2.2). 

In doing so, we shall find it useful to keep in mind some electrostatic analogs of the 
equations we happen to be dealing with. These emerge clearly if we write the current 
conservation equation at the site k as 


jtotal 


4’+/?,. + /S 


( 2 . 6 ) 


Here ^ = 1,2 represents the x and y directions respectively. Furthermore, 

= sinA^^fc, = A^0j and IQ are the total, superconducting, normal and 
external currents, respectively, at node k in the /j*'* direction. Using the TCC condition, 
V • = 0, we have 


= + (2.7) 

On comparing Eq.(2.7) with the Laplace Equation VV = -ptotail^o and recalling 
that p total = Pfree + Pbound, we axrive at the following correspondences 


Tgq- 

III 

(2.8) 

-P(f) = f! 

(2.9) 

III 

(2.10) 


Thus the condition imposed above on the 0,- is nothing other than a choice of 
reference potential or gauge. We, furthermore, note that the JJA thought of as a 
dielectric medium is highly non-linear. Indeed, = sin{- J Ek,,.dt), where = 

Ajk is the component of the electric field at site k. This is to be contrasted with 
P{f) oc E{f) for a linear dielectric. 
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2.2.2 The O(NlnN) procedure 

Going back to Eq.(2.3), can alternatively be inverted by defining the Green’s 
function [see Jackson] 

G{F,n = E f 

Aj;#0 '^k 

where Ag are the eigenvalues and the orthonormalised eigenvectors of Gq^, i.e. 




( 2 . 12 ) 


and f are points on the lattice. More explicitly, the site coordinates (x,, y.) = fi 
referred to an origin located at the lower left hand corner of the lattice, are related to 
the site indices i used earlier by i = (x, — xq) + {yi — yo)Nx + 1 where (xq, yo) denote 
the coordinates of the lower left corner of the array. We note that to get a well-defined 
function, we must explicitly remove from the sum in Eq.(2.11), the zero mode, which 
necessarily exists, since Gq ^ is singular. 

The removal of this mode, however, creates a slight problem. Indeed, if we premul- 
tiply by Gq^ and use Eq.(2.12) we get 


G„-' (f, f")G(r", r) = i; ') 

r " 

r ‘ ■ 

= d[r — r } 


(2.13) 




V V 


^ LxL^ 


(1. 1.--1) 


which in discretised form is 




1 

LjjXjy 


(2.14) 


Now if we multiply Eq.(2.4) from the left by G we get 6, - which reads just 
right, viz. 


I 
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ei = Gijdj (2.15) 

provided the 0,- satisfy = 0. In other words, with an appropriate gauge choice, G 
does allow us to invert Gq^. Indeed, the requirement that it does so fixes the gauge 
uniquely. It is to be noted that Eq.(2.15) involves [d] and not [D], 

The function, (?(f, r'), can be easily evaluated for a periodic lattice. The form of 
the operator case, is given by 

^0 (^5^ ) ^r+ey,r' ^r— et,,r' (2.16) 

It is easily checked that the normalised eigenvectors are of the form and corre- 

LxLy 

spond to the eigenvalues 


= 4 — 2 cos kx — 2 cos ky 


( 2 . 17 ) 


We note that Ag_Q = 0. 

Due to the rectangular periodicity of the lattice the eigenvectors must satisfy 
= '^fc(^+ Ltit) i = x,y. To fulfill this condition, the wavevectors must in 
turn be quantised as k = {2rvK)lLi i = x,y and n = 0, 1, • • • (L,_i). The Green’s 
function of interest is consequently 


(2.18) 


We can now estimate the number of steps required to evaluate G{f,r')d{f'). 
The sum over r' is a fourier transform and can be carried out in NlnN steps to 
produce d{k). Then come the N multiplications \%d{k) followed by the sum over k, 
which has the form of an inverse fourier transform, and requires an additional NlnN 
steps, making a grand total of N{2lnN + 1) multiplications. 

Quite generally then, 9(f) can be evaluated from Eq.(2.15) by (see Fig. 2.1) 

(i) creating an constant x Ly matrix Go(k) = Nj;/X^ where N^ are the normali- 
sation factors due to 
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(ii) e\'aluating d(k) = d[r')'^j^(r'), this being the forward transform (W), 

(iii) computing d'{k) = Go{k)d(k), and 

(iv) taking the 2-D inverse transform W to get 6{t) = d' (k)'^ j;{r) . 

For certain types of transforms and certain values of and Ly, W and W can each 
be performed by matrix decomposition or row-column techniques with a complexity 
0{N In N). We now turn to a discussion of this procedure in the context of finite 
arrays. 


2.2.3 The Imposition of Boundaries 

For a finite array, the general form of (Eq.(2.11)) is modified since now fewer than 
four bonds meet at all sites on the boundary (see Fig. 2.2.3). For non-corner points on 
the left edge of the array for example, the operator has the form —Sr±iy,r' ■ 

However, if ^ satisfies ^(xo,y) — ^^(a:o — l,y) = 0 V y along this edge, one can 
continue using the periodic form of Gq^ because then — ^rdbe^.f' — v) 

automatically reduces to {36^,?! — y). We thus see that the 

finiteness of the array imposes the following boundary conditions on ^ where (xq, yo) 
and {xl, yi,) are the diagonally opposite corners of the rectangular lattice. It is amusing 
to note that the left-hand side of Eqs.(2.19) are the discrete derivatives of ^(x, y) along 
the edges x = xq, y = yo^ x = xl and y = yt^ respectively. Thus the finite case poses 


for us a discretised Neumann boundary value problem. 

^(xo,y) - - l,y) = 0 Vy (2.19a) 

^(x,yo)-$(x,yo-l) = 0 Vx (2.19b) 

y) - + 1, y) = 0 V y where xi, = Xq + - 1 (2.19c) 

'3?(x, yi) - yn + 1) = 0 V X where yz, = yo + Ty - 1 (2.19d) 
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2-D Forward 
Transform 
{W) 


Multiplication 
with G{k) 


2-D Inverse 
Transform 
{W) 



Figure 2.1: The general flowchart of the 0{NhiN) method. Given d{x,y), one per- 
forms a 2-D forward transform and multiplies the resulting d{k^, ky) with Go(A:^, ky) 
before taking a 2-D inverse transform to get d(x^ y). 
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Figure 2.2: An finite Lx x Ly array with current being injected at the left edge and 
extracted from the right edge 



Furthermore, we note that (see Eq.(2.17)) has the symmetries of the square, 
i-e. Xk^^ky — X±k:,,:fky = X.k^^ky- Thus any linear combination of the corresponding 
eigenvectors. 


y) = oi exp i{kxX + kyy) + 02 expi{—kxX + kyy) (2.20) 

+03 exp i{kxX — kyy) + 04 expi{—kxX — kyy) 

is an eigenvector of with the same eigenvalue Ajj. To find the specific linear 
combination, which satisfies the boundary conditions given above, we need merely 
impose each of Eqs.(2.19a-2.19d) in turn. Eq.(2.19a) can be satisfied by choosing 
Xo = 1/2, at = 02, 03 = 04 and the resulting wavefunction is 

T^(a;, y) = 2oi cos kxX exp(ikyy) + 203 cos kxX exp{—ikyy) (2.21) 

Subjecting Eq.(2.21) next to the constraint Eq.(2.19b), we get 


^f(x, y) = 4oi cos kxX coskyy 


( 2 . 22 ) 
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where yo is now required to be 1/2 as well. The fact that xg = yo = 1/2 means that 
the origin of coordinates is chosen on the dual lattice and both coordinates of every 
lattice point are half integers. 

The boundary conditions Eq.(2.19c) and Eq.(2.19d) at the right and upper edges 
of the array can be satisfied by imposing a quantisation condition on k. Indeed, using 
Eq.(2.19c) we have 


-Soi cos kyy 5m{ka;L:c) sin(-;^) = 0 V y 

whereby k^ = UxTc/Lx-, rix — 0, 1, • • •, Tj. — 1. Similarly using Eq.(2.19d) 

UyTK/Ly, nj, = 0,1, • - 1. 

Thus finally, the orthonormalised wavefunctions we are looking for are 


(2.23) 

ky = 






11 

1 - ) cos{kxx) cos{kyy) k^O (2.24) 


2 2 

The resulting Green’s function to be used in Eq.(2.11) is thus[122] 
4_ 

r 


“<*^**') ( 2 . 25 ) 


The action of this Green’s function on a uniform current drive, y) = 

can be rigorously shown to be (see Appendix 2.A ) 


■'X.Xo 


Y^G{x,y, x', y')I{x',y’) 
1/ 



(2.26) 


2.3 Inclusion of Bus-bars 

Using the techniques of the previous section we generalise the procedure to arrays with 
either a single or a double busbar (henceforth referred to as SBB and DBB respectively). 

In the SBB case, the busbar is normally placed along the current extraction edge 
(see Fig. 2.3). The corresponding d and 6 are zero for all time i.e. they do not 
evolve. The drive edge is however open and can be connected to any kind of current 
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Figure 2.3: An array with a single busbar at xl — xq + Lx — 1. The injection edge Xq 
can be driven by any current profile. 


profile. Such an array, with a linear profile, has been used previously to study vortex 
dynamics[78, 79]. On the other hand, an DBB array has shorts along a pair of parallel 
edges (say those at rc = Xq and x = xj, + 1 as in Fig. 2.3). In this case, no current can 
be injected in the x- direction but the y- direction is available to an arbitrary current 
drive. Such an array with electrically connected busbars could be used to simulate e.g. 
the “channel” in the experiments conducted by Van-der-Zant et aZ.[3]. 

We note that by setting the 6 along the busbar to zero and measuring all other 
O's with respect to these preassigned variables, we are unambiguously fixing the gauge. 
Furthermore, since the wavefunction in a rectangular system is separable, i.e. ^{x, y) = 
'^i{x)'^ 2 iy), any change in the boundary conditions along the x- direction, has no 
impact on '^ 2 ( 2 /), which therefore continues to have the form derived in Sec. 2.2.3, viz. 


My) = 


1--4 

2 


cos kyy 


The altered b.c.s however enter crucially into the determination of ’^i(x) and hence 
in to that of ’^'(x, y). 
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2.3.1 Arrays with a Single Bus-bar 

For the SBB case, we have the following b.c.s 

^i(xo) - - 1) = 0 Vy (2.27a) 

'Fi(xi + 1)=0 Vy where xi, = Xo + Lx - 1 (2.27b) 


The second of these results from the fact that the wavefunction must vanish at 
the shorted edge as discussed above. Condition Eq.(2.27a) fixes the form of ^i(x) 
to be ~ cosfcjX where x is a half-odd integer, while Eq.(2.27b) quantises k^: = 

((2^^x + I)’’’)/ (2Lx + 1) Tia, = 0, 1, • • • , Lx — 1. The normalisation constant is fixed by 
the requirement that cos^ fcx(p “ 1) ~ Thus 


%{x) = 


cos kxX 


VLx + 

and hence the Green’s function G for the case of a single busbar is 

~ cos(fej;x) cos{kyy) cos{kxx') cos{kyy') 

4 — 2 cos kx — 2 cos k^ 


(2.28) 


(5(r, f') = 


uE- 


T„(Lx + 1) ^ 


(2.29) 
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The action of G on a uniform input drive can, as in the case 

of Eq.(2.26), be exactly shown to be (see 2.5) 


J^G(x,y, x',y')I{x', y') 



(2.30) 


Although a ID Fast Cosine transform can be used along the y- axis, this cannot 
be done for the x-direction, since no fast Cosine transform of the form 


L-l 

C(m) = ^ C{n) cos 

n=0 


(2n + l)(2m + l)7r 
(41 + 2) 


is known to exist[127]. If we resort to usual matrix multiplication along fia,, the algo- 
rithm becomes 0{2NNa: + 2N\a.Ny-\-N) in complexity and is useful only if Ny » N^- 
A much more efficient approach is to combine the FCT along Cy with cyclic reduction 
along the x-direction. This can be carried out as follows[128]. We begin by writing 
Eq.(2.4) as 






_L a 


(2.31) 


where the operator Ty'^y is defined to be T^^y = V.y-i + ^y'.y+i “ Taking the 

Cosine transform along the y- direction i. e. setting dx,y = ^x,ky cosfcyy, we find 
that 


^ 2(COS ky - 2)dx,ky COS 

ky 

As a consequence the y— cosine transform of satisfies 

^x-l,ky + ■^(^y)^a:,fcs, + = —dx,ky 


(2.32) 


(2.33) 


where we have used X{ky) = 2(cosfcy - 2). Writing Eq.(2.33) thrice with x set equal to 
X - 1, X and X + 1 respectively, multiplying the second of these equations by -A(^y) 
and adding all the three of them together, we get 


Ox-2, ky + ^^^\K)0x,ky + 0x+2,ky = -dx-l,ky + KK)0’XjCy dx+l,ky 


( 1 ) 


(2.34) 
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where = 2-{X(ky))^. We note that the x values occurring in Eq.(2.34) increase 

in steps of 2 as opposed to 1 in Eq.(2.33). The two equations are, however, of the same 
form. We can thus use the reduction procedure repeatedly to get equations with x- 
values increasing in steps of 4, 8, • • • until after m stages (2”* = L^), we end up with a 
single equation for the central line of variables: 

^(m ^^(^ky)9xa+Lxl2,ky = “ ^^o,ky ~ ^xi+l,ky (2.35) 

The and occurring in these equation are defined by the recursive 

relations: 

X^^\ky) = 2-(x^^-^\ky)y (2.36) 


Jp) _ 

“r.fcy ~ “a:-2P-i,fcy 




(2.37) 


Thus, if and OxL+i,ky are known, then 9xo+Lx/2,ky can be immediately deduced. 
The knowledge of ^xo+L*/ 2 ,iy leads to that of dxo+Lx/i,ky and 9xo+3Lx/4,ky which in turn 
determine 9xa+qLxls,ky (g = 1, 3, 5, 7) etc. As for the starting values, ^xr+i.fcy = 0 V fcy, 
since this is the y-cosine transform of 9x^y for x — xj^ + 1, i.e. taken along the busbar. 
9xo,ky must however be determined by explicit multiplication. Fortunately, since we 
know the relevant Green’s function explicitly, we can carry this out in 0{N) steps, by 
taking the cosine transform of Eq.(2.15): 


'XQ ,ky 


Ly{L^ + I) 


EE 


(l - 


4 — 2 cos kj. — 2 cos k. 


cosky;XQCosky;x'dy:>^ky (2.38) 


The number of steps required to evaluate all the higher-order divergences d^^\ , p = 
1, • • • m— 1 used by the method can be performed as follows. For a given Nx the number 
of higher order divergences that need to be calculated (starting firom level 0 which 
correspond to ordinary divergence) is given by / = ln 2 Nx - 1. At the (/ - m)*^ level 
there 2”"+^ - 1 divergence vectors which need to be calculated. Each such calculation 
requires Ny multiplications. Thus the total number of multiplications is 
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ii (2’”''' - 1) 


m=0 


N,. 


L 2 


- In2 iV, 


(2.39) 


With these divergences in hand, the cyclic determination of 6x,ky tas a complexity 
0{N). Adding to all this, the number of multiplications involved in carrying out the 
forward and backward FCTs in the y-direction, we conclude that the entire procedure 
can be carried out in 0(2 A In Ny + 2N + N/2 — Ny Ina Nx) steps. The restriction of Lx 
to integral powers of 2 can be removed as shown in Ref. [129]. The whole procedure is 
described in Fig. 2.5. 


2.3.2 The Case of a Double Bus-bcir 

Turning now to the DBB case, we note that the wave-function, ^'i(x), must now satisfy 
the following b.c.s 

^a(xo) = 0 yy (2-40a) 

^i(a:L + l) = 0 Vy (2.40b) 


Condition Eq.(2.40a) fixes the form of ^i{x) to be ~ sinkxX where x is an integer 
while condition Eq.(2.40b) restricts kx to the values {nx-K)/{Lx + 1), = 1, • --Lx- 

The normalisation constant. A, is fixed by A^Y^^^siiP kxP = 1. The resulting or- 


thonormalised wavefunction is then 


^i(x) = 



sin kxX 


(2.41) 


and hence G for an array with a double bus bar is 

4 n -Ui. .)sm(kxx) cosik^y)sm{kxx')cos{kyy’) . 

G{f,r') = 4 - 2cosfci - 2cosfcj, 

This G coincides with the Go for this case because the zero-mode is absent. In this 
case, the transformations along both directions are amenable to fast, i.e, O(mnN), 
procedures provided + 1 end Ly are both integral powers of 2. 
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dx^y 



dx,y 


Figure 2.5: The general flowchart of the procedure described in Sec.2.3.1. Given 
one performs a 1-D forward transform to get d^^ky, calculates the higher order diver- 
gences (d^^ky) tridiagonalisation to get Ox,ky Finally, a 1-D inverse transform 

provides 
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For the case of periodic boundary conditions in one direction (say y), one can use 
the fast algorithms by simply noting that ^^ 2 (y) ~ expife^y with ky = 2 n 7 r/Ly. In this 
case, one would use a fourier as opposed to a cosine transform in the y direction, for 
the SBB case. For the DBB case, cyclic reduction can be used in the x direction with 
a fourier transform along the y or alternatively in the y direction with a sine transform 
along X. Arrays with such boundary conditions have been used by m any authors in 
their numerical studies[57, 74, 77, 109] 

2.4 Defects 

We move on next onto arrays containing defects in the form of broken bonds. As 
mentioned in Sec. 2.1, arrays of this sort have recently attracted the attention of several 
groups. Xia and Leath showed that in contrast to networks of resistors[l30, 131], where 
breakdown emanates from the most critical defect, JJAs with missing bonds do not 
become resistive, the moment any one junction turns critical. Cohn et a/. [104] showed 
that defects exert an additional pinning force on a vortex placed in the system. In 
Chap. 4 we demonstrate the existence of multiple vortex sectors in the steady state 
configurations of the system. This is done by running on large arrays using a modified 
version of the above algorithm which we now describe in detail. 

The evolution equation 9 for an array containing a single missing bond between 
the sites k = (xq, yo) and (k + l) = (xq + 1, yo) can be written in matrix notation (see 
Eq.(2.2)) as 

(Go'+/^)[0] = [d] (2.43) 

where Gq^ is the discrete Laplacian with free, periodic, SBB or DBB boundary condi- 
tions and h is given by 

hfcjfc ~ “ 1 

hk,k+i = hk+i,k = 1 

hi,j= 0 i,j^k,k + l 


(2.44) 
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Multiplying Eq.(2.43) with G from the left we get 

I 

(I + A)[e]=^G[d]=^ (2.45) 

where A = Gh and I is the identity matrix. The process of determining G[d] = ^, as 
already been outlined and can be executed in N In N steps. One then needs to evaluate 
(1 + eflBciently. To this end, we note that A has the form 

Aij = (2.46) 

where the numbers Xi are got by explicitly by multiplying Go and h. From Eq.(2.46) 
it follows that 


A'' = A(Ae)«-^ (2.47) 

/ 

where Ae = Xk — Xk+i- In fact, Ae is the only non-zero eigen value of the matrix A 
and for the case of SBB discussed in Sec. 2.2.3 it is given by 

Ae = -Azr- S i ^ 

JL/ JL/ fj 4 2i cos rh^ “* 2i cos fc*# 

^ ^ k^O ^ 

where ki = {niTr)ILi^ n, = 0 . . .Li — 1 and i — x,y. From Eq.(2.47) it follows that 

(I + A)-^ =I-A + A^-A^ + --‘ = I- — ^ (2.49) 

^ 1 -I- Ae 

Since the constant matrix R= — A/(l-f Ae) contains only two non-zero columns, 

which are moreover negatives of each other, [^] can be determined from ^ in just N 

multiplicative steps. 

The extension of the above procedure to the case of a linear defect[77] consisting of 
n broken bonds is straight forward as each additional defect introduces into h a 2 x 2 
diagonal block of the form given in Eq.(2.44). In this case, the series (Eq.(2.49)) cannot 
be summed analytically due to the presence of cross terms. An analytic summation is, 
however, unnecessary. It is sufficient to note that 


Ay = — A*2+i) •! - Afcn+i) (2.50) 
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for bonds missing between fc,- and ki^i (i = 1,2 — n), and that consequently A'* has 
the same form as above. The numbers Xf’s occurring in A'^ are of course 9— dependent. 
It follows from this that the series can once again be written in the form 

(J + A)"^ = J + (2.51) 

where is a constant matrix of the form given in Eq.(2.50). Clearly the action of R 
on a vector of length N can be determined in multiplicative steps. The overall 
complexity of the algorithm in this case is thus 0(N In N + nA^). Since the num- 
ber of broken bonds required to observe various breakdown phenomena[77] is much 
smaller than N, this results in a very large saving in the time required to perform the 
computations. 

It is interesting to note that the Green’s function for the broken bond system 
satisfies the Dyson’s equation[132] 

G** = G -f G(-h)G^ (2.52) 

The matrix (—h) is thus the analog of the potential due to an impurity in an 
otherwise perfect lattice and the case of the single defect is the analog of the S function 
impurity at a given site in the system. 

The procedure developed above is for the situation where no pair of broken bonds 
has any site in common. This ensures that the block matrices introduced in k are 

always 2x2 rather than m x m, m > 2 The general case can be dealt with in a 

« 

similar manner albeit with an increase in complexity. Lastly, bonds can be added in 
the interior of the network rather than being eliminated since this does not violate the 
constraint T2i d,- = 0. 

2.5 Summary And Discussion 

To summarize, we have extended the Eikman-Himbergen algorithm to a number of 
practically occurring situations. In particular, we have extended it to the case of a 
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bus-bar placed along one edge of a rectangular array as also to that of electrically 
connected bus-bars placed along two parallel edges. Finally we have shown that our 
simulation can be carried out, with only a marginal increase in complexity, for bonds 
eliminated from or added to the array. In all cases, the open edges can be connected 
to current sources having any profile whatsoever. 

The algorithms for bus-bars as also for perfect arrays subject to magnetic fields, 
noise, and/or current drives, all involve Green’s functions specific to the lattice in 
question. The Green’s function of interest can straight-forwardly be constructed once 
the eigen-basis of the corresponding connectivity matrix or discrete Laplacian has been 
determined. We have shown that this eigen-basis consists precisely of those linear 
combinations of eigenfunctions of the connectivity matrix for a periodic lattice, which 
satisfy an appropriate set of boundary conditions. 

The addition and removal of bonds introduces additional boundaries often in the 
interior of the array. As a result our earlier procedure becomes invalid. An island of 
defects created at the center of the array, e.g., defines a lattice version of the sinai 
billiard problem, which classically has only chaotic solutions. However if the number 
of bonds eliminated is small. Green’s function techniques which have been extensively 
used to study disordered solids can be very elfectively applied. We have shown that this 
techniques 3delds an 0{N In N + nN) algorithm. The extensions of such algorithms to 
the case of 3D arrays is easy since the corresponding eigen-vectors are specified to be 
either exp(ifc ■ F) or cos(kxx) depending on whether the boundaries are periodic or free 
boundaries respectively. The eigen-values are now changed to 6 — 2coskx — 2 cos ky — 
2coskx. In this case, one has to obviously use 3-D transforms. 

It is also worth noting that many of the expressions we have derived can be written 
in terms of continuous functions, which are, of course, to be evaluated or sampled only 
at points belonging to the array. Now if we turn this observation around and think of 
the discrete lattice as being produced by the sampling of the continuous 2-D waveform, 
we make contact with a long-standing problem in electrical engineering, viz. that of the 
digitisation and subsequent recovery of band-width limited analog wave-forms. Usually 
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such waveforms are processed as rectangular spaced arrays i. e. they are periodically 
sampled along the two orthogonally independent variables. However, one can also use 
hexagonal sampling[133]. Once discretised, the waveforms can be processed as arrays 
of numbers x(ni.n 2 ) by the computer where ( 711 , 7 x 2 ) is a discrete point on which the 
observable x has some value. The periodicity of the observable x( 7 Xi, 7 i 2 ) is important 
as it specifies the eigenfunctions associated with the system. The matrix has 
precisely the information about the chosen sampling raster and the periodicity of the 
lattice. The raster defines the total number of nearest neighbours of a general site 
( 711 , 7 x 2 ) while the periodicity (or finiteness) of the lattice decides the connectivity of 
the sites at the boundary of the network. The eigenfunctions used in both cases are 
the same. Indeed, the analogy goes deeper and we can think of the JJA itself as 
a latticized version of the continuum case (superconducting islands embedded in a 
normal background much as the Abrikosov lattice consists of normal regions embedded 
in a superconducting background). Finally, since hexagonally sampled waveforms can 
be recovered in 0{NlnN) steps, this analogy opens up the very interesting possibility 
of working out 0{N In N) algorithm for a triangular JJA (see Chap, 3). 


Appendix 2. A 


In this appendix we show that the effect of a constant current applied at the edges of 
a finite array. The sum as in Eq.(2.26) can be written as 


X{x, y) = / E G{f, x' = 1/2, y') - G(f, x' = L, - 1/2, y') 

.y' 

4/ __(l-l/2<5fc,,o-l/24,,o) , 

= > > ; T ; — - cos k^x cos kyy 

L^rLy ^ ^ 4 — 2 cos kx — 2 cos ky 


(2.A.1) 


cos kx / 2 — cos kx{Lx — -) cos kyy' 

Noting that Y^y> cos kyy' = LySk^fi and fca: ^ 0 (since ky = 0) this gives 

- / cos sin^ cos ^ 

X{x) = y- E r- 2 ir 

k^^O ■# 


(2.A.2) 
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independent of y. Since x, y and y' are half-odd integers and k = (n,' 7 r)/L.-, i = x, y 
we can rewrite the above expression as 


V 21 '7r7l(77l + 1 /2) 

X(m) X{n)cos—^ LI 

■ Lj't 


n=l 


(2.A.3) 


where 0 < m < {L^ — 1) and 


X{n) 


Siscos ^ 


•sm 


2 nTT 

2L^ 


Now we claim that 


X(m) = I (^2^ 

To prove this we define A(0) = 0, then 


- m = / 


f Lx 


\ 2 


■ X 


(2.A.4) 


(2.A.5) 


X(m) = cos 


X(n) 


where ko — l/\/2 and — 1 for n 0. This is exactly the forward cosine transform 
of A(n)[127]. The inverse cosine transform is 


xr/ \ 7 -{>/ \ (2771 -1- liriTT 

A (77) = ^ X{m) cos ^ ’ 

m=0 


2L, 


which for tt = 0 gives 


1 1 




(Tx- 1)1, L,(i, - 1) 

2 2 


= 0 


as required. 

For 77 ^ 0, we get ^( 77 ) = m cos since the first term yields 0 

when 77 0. Expanding the cosine term and using the identities[134] 

mkir sin Aitt T, cos 


E 


777 sm 


4 sin 


2 kit 




m=0 ^x ^ 

mk-K 7, sin i_cos%is 


E 

m=0 


mcos 




4sin^fe 


(2.A.6) 



CHAPTER 2. FAST ALGORITHMS FOR SQUARE ARRAYS 


41 


This yields the required result for X(n) (Eq.(2.A.4)). 

This appendix shows the effect of a constant current applied to a SBB configuration. 
To evaluate the sum as in Eq.(2.30) we once again use the above procedure whereby 


X(x,!/) = /^G(F,i'=1/2.s') 

^ (l ~ 

(Lx + l/2)Lj, “ 4 - 2 cos fca; - 2 cos ky 

cos kxX cos ky-y cos kx/2 cos kyy' 


(2.B.1) 

(2.B.2) 


where ky = (nyTT)fLy and k^ = ((2na; + l)7r) /{2Lx + 1). Noting once again using the 
fact the the y' summation yields LySk^.o and that x, y are both half-odd integers we 
can rewrite the above expression as 


X(m) = 


21 

Lx + 1/2 


Lx — 1 

X(n)cos 

n=0 


(2n + l)(2m + l)7r 
ALx + 2 


where 


X{n) 



(2n + l)7r \ 

ALx + 2 j 



(2n + l)7r '\ 
ALx + 2 J 


We claim that 


X(m) = I(Lx -n)= I(Lx -^ + \) 


(2.B.3) 


(2.B.4) 


(2.B.5) 


To prove this one needs to first determine the inverse transform of Eq.(2.B.3). Although 
such a transform is rarely used, one can show that the inverse transform is 


XT/ ^ <rf \ (2n + l)(2Tn + l)7r 

X(n)= E X(m)cos 


m=0 


(2.B.6) 


Now using Eq.(2.B.5) in Eq-(2.B.6) cind the identities (Eq-(2.A.6)), it can be shown 
that indeed one gets Eq.(2.B.4). 



Chapter 3 


Fast Algorithms For Triangular 
Arrays 


3.1 Introduction 


Various experimental[32] and numerical studies[30, 37, 74, 76] studies on JJAs indi- 
cate that vortices, which are natural excitations of the system, play an important role 
in determining both the equilibrium and dynamical properties of these arrays. The 
behaviour of vortices, in turn, depends in large part on the vortex pinning potential, 
which arises from the discreteness of the lattice. Many features of vortex behaviour 
are thus better understood by varying the pinning potential. This can be achieved 
by switching to a triangular lattice, which has a ~ 5 times smaller pinning potential 
than a square array [61]. Triangular lattices are hence better suited for the study of the 
dynamics of an ideal 2D gas of vortices[4]. They also provide a better chance for the 
observation of ballistic vortex motion, i.e. the motion of a “massive” vortex through a 
force-free region due to its own inertia, which has been predicted theoretically [135, 55, 
53, 136] and has to some extent been investigated experimentally[3] and numerically[56, 
57, 59, 137]. The geometry of the hexagon, which occurs naturally in triangular arrays, 
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allows us to conveniently study the Aharanov-Casher effect wherein vortices, as quan- 
tum objects, can interfere constructively or destructively after following two different 
paths enclosing a charge placed at the center of the hexagon[119]. The ground state 
configuration of triangular arrays under the application of an external magnetic field 
has unique features of its own. For example domain wall superlattices[138] occur at 
characteristic values of /, the quantum of magnetic flux threading a plaquette, and 
there is an accidental degeneracy of the ground state for / = 1/3 and 1/4[139]. In 
fact, there are indications that the case of / = 1/4 belongs to the same universality 
class as the fully frustrated / = 1/2 case of the square lattice[94, 140]. It has also been 
observed that the presence or absence of fractional Shapiro steps depends on the rela- 
tive orientation of the lattice and the external current drive[l41]. Lastly, we note that 
a triangular configuration of three capacitive [142] and resistive[143] junctions offers 
the simplest geometry capable of sustaining a vortex while being amenable to explicit 
theoretical analysis. 

As indicated in Chap. 2, a number of authors have evolved fast algorithms for 
square lattices with a variety of boundary conditions. No fast algorithms have however 
been discovered so far for the case of triangular arrays making numerical work with the 
latter prohibitively time consuming. We, accordingly, take up next the task of filling 
this lacuna. 


3*2 Notation 

The algorithms which we develop in this chapter can very conveniently be discussed 
in terms of the bracket formalism of Dirac[144, 145]. A given state of the system 
corresponds to a specific configuration of the ^’s and can be denoted by \6), In the 
position basis {|x) ~ |xi,X 2 )}, the components of \d) are obtained by taking the usual 
projections {xi^X 2 \d) = 0{xi^X2)- We note that the coordinates (xi,X 2 ) of the various 
lattice sites are related to the labels i used in Eq.(2.4) by (x— xo) + (j/-~t/o)Aa: + l where 
(xo, t/o) are the coordinates the bottom-left corner of the array. We could alternatively 
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consider the components 6 {ki,k 2 ) = in the momentum basis These are 

related to 9 {xi^X 2 ) by the unitary transformation matrix (fc|f): 

^2) = (3-1) 

X 

where the summation is to be carried out over all the lattice points. The validity of 
this equation as also the unitarity of (^js) rests on the completeness of the two bases 
in the N- dimensional space of interest: 

j = y:ix>(i|=y;ii)(*i (3.2) 

X k 

where I is the N dimensional Identity matrix. Both bases are further assumed to be 
orthonormal, i. e. 


{x\x') and (3.3) 

We note that the delta functions occurring in Eq.(3.3) are the kronecker (as opposed 
to dirac) deltas. Furthermore, Eq.(3.3) indicates that in the position basis, the vector 
|ic) is zero at all points other than x = (xi,X 2 ) where it is 1. 

The 6 i and the divergence d,- similarly give rise to the vectors |^) and |d). Eq.(2.15) 
is then 

<?„-'!«) = l<i) (3.4) 


written out in the position basis. 

The underpinnings of the fast algorithms can now be made transparent by writing 
Eq.(2.15) as follows: 

{x\ 6 ) = {x\G\d) (3.5) 

k k’ 


We note that the r. h. s. of Eq.(3.5) involves 
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(i) Evaluating This is the forward transform W and can 

be executed in certain cases with complexity O(JVlnjV) 

(ii) Multiplying {k’\d') by {^G\k') i. e. the operator G expressed in the k basis, to 
get {k\df). If G happens to be diagonal in this basis, the multiplication can be 
carried out in N steps. 

(iii) Taking the backward transform W i. e. Here too, the complexity, 

for some cases, is O(NljiN). 


It follows that we have a fast algorithm for a given lattice as soon as we have dis- 
covered a basis {\k)} which (i) is related to the position basis {jx)} by a fast transform 
of some sort and (ii) completely or almost completely diagonalises the corresponding 
G or Gq ^ . 

For a periodic Li x L 2 square lattice (where Li and L 2 refer to the lengths in the 
X and y directions respectively) 

(x^|C?Q^|x) = (45x,i' ^x,x'+ei ^x,x'—ei ea ) (^•®) 


and the basis defined by 


(xl*) 


1 


:exp 


^ik ■ xj 


(3.7) 


y/LiL2 

with ki = 2 nj 7 r/E,-, does the job. The transform W and W in this case are the usual 
direct and inverse fast fourier transforms respectively, and G is clearly diagonal: 

1 1 


{k'\G\k) = 


ik'\k) = 


^kx fky 


4 — 2 cos kx — 2 cos ky 
For a finite square lattice, the Go ^ matrix, on the other hand, is 


(3.8) 


(x'|(?Q ^ |x) = ei )(1 d" (^5,x' ^x,x'+ei)(l ^xi .xo+ii — 1 ) 

+(^x,x' ~ ^x,i'-e2)(l ~ ^x2,yo) d" ( 5 x,x' “ ^x,x'+e2)(l “ l) 

where (xq, yo) are the coordinates of the bottom-left edge of the array. The k- basis 
is now given by 





{x\k) = iVg cos AiiXi cos ^ 2 X 2 (3.10) 

( 

where ki = rnrfLi, 1/2 < Xi < Li — 1/2; i = 1,2 and iVj are the normalisation con- 
stants. We note that the set of [fc) vectors spans the vector space of states/configurations 
and moreover diagonalises G. The transforms W and W are now the direct and inverse 
fast discrete cosine transforms. 

It is to be noted that the {\k)} basis is different for different lattices. Strictly 
speaking, we should introduce a new symbol each time to indicate this. However, 
that would make the notation cumbersome. We, accordingly, use the symbol {[fc)} 
uniformly and let the context and the explicitly stated form of (x|A:) decide its precise 
meaning. 

3.3 Periodic Triangular Lattice : Sampling Strate- 
gies 

In contrast to square lattices which have four nearest neighbours, triangular lattices 
have six neaxest neighbours. The connectivity matrix is therefore different in this 
case. Moreover, the bonds are at a relative angle of 60° to each other. Such configu- 
rations are found in crystallographic studies and in various radio physics applications 
(e.g. radar). We, therefore, digress for a moment to study how such lattices are con- 
structed especially in the context of Digital Signal Processing. In the process we learn 
some important techniques applicable to the JJA context. 

It is well known that both square as well as triangular lattices can be created by 
sampling a 2-D continuous waveform can be sampled at certain discrete points[146]. 
The value of any variable, say $, is measured at these discretised points. Intuitively, 
rectangular sampling i.e. periodic sampling along the two orthogonally independent 
directions is the most common. Such a sampling produces, in general, a rectangular 
lattice and hence a discrete $ field. 


(n|$) = $(ni, 712} = $a(nia, n2b) where Ui, n2 = 0, 1, 2, • • • 


(3.11) 
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and is the corresponding continuous (analog) waveform. One of the many assets of 
using such sampling is its separability whereby the sampling along each of the axis is 
independent 

(n|$) = (ni|$i)(n 2 |$ 2 ) = (3.12) 

Here a and b are the units in which sampling is done along the two directions. 

There exists other sampling strategies e. g. hexagonal sampling which essentially 
produces a triangular lattice (see Fig. 3.1). Mathematically, this sampling is described 
as 

{ff|4) = $(ni. n,) = n,b) (3.13) 

whereby each odd row is shifted half a unit from the even rows in the rij direction. 
Unlike rectangular sampling, this is not separable in the variables ni and n 2 - However, 
as shown by Merserau[133], such a sampling has the added ad\’antage in that it requires 
fewer samples to describe the analog waveform restricted to a particular area. 

Periodicity of the observable $(7Zi, 712 ) is another important concept and finiteness 
(as in finite JJAs) essentially disrupts the periodicity. In the most general case, the 
observable $ is said to be periodic if it obeys 

$(n) = $(n + iVr) (3.14) 

where iV is a 2 x 2 matrix and r is any 2D integer vector. Here too, the most obvious 
choice is rectangular periodicity whereby 

$(ni, 712 ) = $(tii + Li, 7 x 2 ) = $(ni, 712 + L 2 ) (3.15) 

However, for some lattices, (e.g. triangular lattices) one can have in addition to rect- 
angular periodicity, hexagonal periodicity which is defined as 

$(7x1, 7X2) = $(ni + 2 Z>i -b L2, 7x2) ( 3 - 16 ) 

= $(7X1 -b Ti + L2, 7X2 + L2) 

= $( 7 X 1 "b Lx, 7 X 2 -b 2 T 2 ) 


I 
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Figure 3.1: A sampling strategy for triangular arrays producing tiles 


Such periodicity, when imposed on an infinite triangular lattice, produce a tiling 
(see Fig. 3.1) of the latter by the array in question (this has been referred to as the 
“fundamental period” in Ref. [133]). This tiling is quite distinct from the one which 
follows from the conditions (3.15). If corresponding points on different tiles are now 
identified with one another, one immediately obtains the connectivity of points on the 
boundaxy of each tile. For example, in the rectangularly periodic lattice of Fig. 3.2, the 
points carrying the same labels on the top and bottom rows (left and right columns) 
are treated as identical. In other words, the lattice is folded like a torus. The folding 
of the hexagonally periodic lattice of Fig. 3.3, on the other hand, is quite unfamiliar 
but follows straight-forwardly from the tiling displayed in Fig. 3.3. The points to 
be identified here have, as before, been given identical labels. We should mention 
that the above identification can be physically be realised through the use of electrical 
shorts. For a given value of Li and L 2 the number of samples present in a tile is 
given by det N. In the case of a diagonal iV, as in the rectangularly periodic case. 
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detiV - L1L2 while for the hexagonally periodic case it is L2(2Li + i,). For further 
details on hexagonal periodicity, we refer the reader to Ref. [133]. Apart from specifying 
the connectivity of an arbitrary site, the underlying periodicity limits the choice of 
eigenfunctions {xlk'j needed to diagonalize the corresponding laplacian Gq^. With a 
knowledge of the sampling strategy and the periodicity of the variable, the important 
task for electrical engineers is the complete recovery of a bandw'idth limited analog 
waveform from the given samples, $ (called the sampling theorem). This is done by 
taking various transforms e. g. Fourier Transform in the case of rectangularly periodic 
samples and Hexagonal Fourier Transforms in the case of hexagonally periodic samples. 
When Z-,’s are restricted to 2^, the corresponding fast transforms are used. The 
matrix in the case of JJAs has precisely the information about the chosen sampling 
and the periodicity (or the finiteness) of the lattice. In fact, the eigen-functions of Gq ^ 
and the transformations used to process such digitised signals are exactly the same. 
It is, thus, natural to believe that JJAs are like the discrete versions of a continuum 
case. Sampling strategies are thus inherently connected to the development of fast 
algorithms in such arrays. 

The above discussion shows that a triangular lattice of Josephson junctions ad- 
mits both rectangular and hexagonal periodicity. Accordingly, we discuss each type of 
periodicity and develop fast algorithms for each in turn. 

3.3.1 Fast Algorithms for Rectangularly Periodic Lattice 

A convenient choice of coordinate system for a rectangularly periodic triangular array 
is shown in Fig. 3.2. The eigenfunctions in this case can be taken to be oc exp{ik • x) 
where k = kiii -f ^1262 and r = xiii 4- 1262 with e-i • 62 = —1/2. The action of Gq ^ on 
these eigenfunctions is best described in terms of ki = ki + fc2/2 and K2 = ^2 + fci/2. 
Indeed, 

(x|Go ^ |x') = — ^x,i'±ei ~ ^S,x'±e2 ~ 5 '±(62-62)) (^• 1 '^) 


whence 






Figure 3.2: A triangular lattice with Li = L 2 = 3 showing rectangular periodicity. 
Equivalent points on the edges are marked with same labels. 

(x|<?o^l^) (3-18) 

x' 

= [6 — 2 cos «i — 2 cos K 2 — 2 cos(/ti — /C 2 )] exp + K 2 X 2 ) 

It is noteworthy that Xk.x,k 2 has all the 12 symmetries of the hexagonal group: 


Xki,K2 ^—Ki,—K2 '^±/C1 ,±(«1 — KJ ) ~''2 ),iK2 (3.19) 

— X^K2,±{ki—K2) ~ -^±/C2i±«l ~ '^±(k1-K2),±Ki, 

Due to rectangular periodicity (Eq.(3.15)), the wavevectors have to satisfy the quan- 
tisation conditions «;,• = (2n,-7r)/L,-, i = 1, 2 with n,- = 0, 1, • • - L,- - 1 and the Green’s 
function for the rectangularly periodic triangular lattice reads, accordingly, 

( 5 (f x') = ^ Y exp i [ftfyxi - x[) + K 2 {x 2 - xQ] 

L 1 L 2 ^ 6 — 2 cos «! — 2 cos K 2 — 2 cos(/ci — K 2 ) 


(3.20) 
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Figure 3.3: A triangular lattice with L = 3 showing hexagonal periodicity. Equivalent 
points on the edges are marked with same labels. 

This is the same as the Green’s function of a periodic square lattice with only the 
eigenvalues changed. Thus a knowledge of the 0(iVln N) algorithm for the square case 
immediately provides a fast algorithm for this configuration, as well. 

3.3.2 Fast Algorithms for Hexagonally Periodic Lattice 

The eigenfunctions of for this case can, once again, be written as exp(ik • x), 
provided we set k = kiii + ^ 2^2 s-iid x = X 2 ^+Xiei with e,- • Cj = 6ij and | “■ 2^2 

(see Fig. 3.3). In terms of ki = ^ki, K 2 = |^ 2 > these functions, on being normalized, 

take the form 

(xlit) = , ^ — expi{Ki(2xi - X 2 ) + K 2 X 2 } (3-21) 

yjL2{2Lr + L 2 ) 

It is easily verified that the action of Gq^ on {x\k) yields 

A .. 

expi{/Ci( 2 xi — X2) + 1^2X2)} 

= iV [6 - 2 cos(ki + k; 2 ) - 2 cos(ki - Kj) - 2 cos(2/Ci)] 


central LiBRARIf 

I. l T.. ICANPtfll 
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expi{Ki{2xx — X2) + ^,2^2)} 

I 

= Akj ^ exp ii{«i ( 2 a;i - X2) + K2X2)} 

where N = [L 2 ( 2 Li + is the appropriate normalization constant The eigen- 

values Xki,k 2 are 12-fold degenerate and the corresponding eigenfunctions are easily 
checked to be 


exp?‘{±«;i(2xi - X 2 ) -t- K 2 X 2 }, expi{±«i(xi +X 2 ) + ^(xi — X 2 )}, 

exp«{±«i(xi — 2x2) -f K2X2} 


(3.23) 


and their complex conjugates. 

All these wavefunctions can be made to respect the hexagonal periodicity of the 
lattice by simply quantising The quantised values of the momenta, for the case 
L\ — 1^2 “ are 


. , TT , n27r 

Ki = (2ni — n 2 )— and K 2 = -y- 
oIj Jj 


(3.24) 


where ni and n 2 are integers, 0 < ni < ZL, 0 < n 2 < T. 

With this in hand, we can write down the Green’s function for the hexagonally 
periodic triangular lattice as 


, 3L-1,L-1 .. 

G(xi,x2lx;,x'2) = — 53 

ni=0,n2=0 

ni,n2^0 


exp 


[^(2x1 - a;2)(2ni - 712 ) -f jX2n2 


(3.25) 


expz 


L3^(2x; - X2)(2ni - na) + ^x^na 


with 


V,n 2 = ^6 - 2 cos(2ni “ ^ 2 )^ - 2 cos( 2 n 2 ~ - 2 cos(ni + na) 


27r\ 

zl) 


It has been shown by Mersereau that one can move back and forth between the x- 
basis and the k- basis defined above in 0{N In N) steps by the forward and backward 
Hexagonally Discrete Fourier Transforms (HDFT) when L = 2"'. This is done by 
noting that the forward transform of d(xi,xa) is a sum 

3L-1 L-1 p ^ TT 

d(ni, na) = 5^ 5D ^(^ 1 -^ 2 ) exp -i yjr(2xi - X 2 )( 2 ni - na) + yXana 

XI =0x2=0 Lo-6 L 
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can be divided into four smaller sums 

d(ni,n2) = Soo(Xi,X2) + See(Xi,X2) + Soe(Xi,X2) + S^o(Xi,X2) (3.26) 

Here the symbols S'oo(xi,X2) etc. signify that both xi and X2 go over odd values (in 
Soe{xi, X2'), Xi goes over odd values while X2 goes over even values). Using the notation 
d(ni,n2) = HDFTn [^(2^1)2:2)] one then immediately gets 

5ee(ni,n2) = HDFTN/2d(2ri, 2 x 2 ) (3.27) 

Seo{ni, n2) = exp - 2n2)HDFTN/2d(2ri, 2x2 + 1 ) 

Stt 

Soe{ni, n2) = exp -z— (2ni - n2)HDFTN/2d(2xi + 1, 2x2) 

27r 

Soo(ni, 712) = exp + n2)HDFTN/2d(2xi + 1,2x2 + 1) 

Knowing the outputs S,, i = (ee), (00), (eo), (oe), one can immediately get d(ni, 912) ^ 

d(ni, 722) = Soo (3.28) 

Q r 

4 ni + — . "2) = s„ + 

d(n, + I, n, + ^ ) = S„ + 

c r r 

<i("i + Y’ "’ + f 


where TUjJ* = exp[— (27rm)/A:]. The calculation is summarized in the flowchart of 
Fig. 3.4. We mention in passing that (xl<i)’s being real leads to a further saving in 
computational time over what has been described by Mersereau. 


3.4 Finite Triangular Lattice 

We now show how boundaries can be introduced without losing the asymptotically 
linear nature of our algorithms. The finite array can be driven by currents injected or 
extracted through the open edges. To keep things as simple as possible^ we discuss a 


I 
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Figure 3.5: A triangular array with = Xj = 3 which is finite in the ii direction. 

lattice which is periodic in the y- direction and open in the x- direction (see Fig. 3.5). 
Square arrays with such boundary conditions have been repeatedly used by many 
authors in numerical studies[77, 109]. 

We can conveniently study the operator for this configuration by separating it 
into two parts, namely, 


Go^=P-B (3.29) 

where P is the rectangularly periodic operator for the triangular lattice given by 
Eq.(3.17) and B is the boundary operator defined by 

^xi,Xo+il-l[2^®-®' “ ^®i+l^'^X2,xJ^Xi+l^J^X2-l,xi] 


(3.30) 
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One can easily check that B deletes all the bonds required to introduce boundaries 
along X = Xo and x = Xq + Z-i — 1. 

We examine next the k basis defined by 


(x|fc) a cosKiXi cosK2(2a:2 + Xi) (3.31) 

where k - + A: 2 e 2 , and K 2 = 1^2, as in Sec.(3.3.2). However, in this 

case X = X 2 C + Xiei with e, • ej = and ^ + |e 2 - 

To motivate this choice, we note that the wavefunctions on the r. h. s. of this equa- 
tion can be constructed as linear combinations of the degenerate hexagonally periodic 
wavefunctions listed in Sec.(3.3.2) and can further be made to satisfy at least in part, 
the appropriate boundary conditions for our configuration (see below). The fe-basis is, 
moreover, accessible from {[x)} through a fast transform. 

The normalisations on the eigenfunctions are easily determined to be 




L1L2 

L1L2 

2 

4 


if = 

if ki = = 0; i = 1 or 2 

if fc,- = ^ 0; i = 1, 2 


(3.32) 


We shall, however, work with the unnormalized eigenfunctions and normalize them at 
the appropriate junctures. 

For these functions to be periodic under X 2 — ^ X 2 -I- L 2 required by the geometry 
of our configuration, we must have K 2 = (m27r)/T2, ^2 = 0, 1, • • • L 2 - 1 but can choose 
X 2 arbitrarily. We shall for convenience work with X 2 = 0, 1, • • •L 2 — 1. 

As before, we can show that the operator P is diagonal in the k basis. Indeed, 


{k\P\k') = 5^(fclx){x|P|:P)(x'|fc') 

= ^(^| x ) - 5xx±1,x[^x2±1,x'^ - 


x',x 


cos k[x[ cos /C2(2X2 + Xj) 


= X^(^|x)(6 -4cos«i COSAC2 - 2cosk2)(x|P) 

X 

) ~ ^Kl,K2^k,k' 


(3.33) 



CHAPTER S. FAST ALGORITHMS FOR TBJANGULAR ARRAYS 


57 


On the other hand, B is only block diagonal, at best. Indeed, {x\B\k') can be 
brought close to being proportional to {x\k') but not quite. More specifically, by 
choosing Xq = 1/2 and k[ = Trm^lLi, m[ = 0, 1, • • - Li - 1, we get 


(x|5|A:^) — 2(l C 0 S/C 2 ) 


xi ,1/1 -“1/2 


{x\9) 


(3.34) 


From this equation it immediately follows that 




(3.35) 


= 2(1 — COS/C 2 ) 


L 2 


/Cl /c( /Cj(2Xfi — 1) /c((2Z'i — 1) 

cos — cos cos — cos 

2 2 2 2 


= 2(1 — COS /Cj) [1 + (— l)'”^'^'"il cos cos 


2Li 2Li 




where = Trmi/Li and k[ = ■nm^/Li. We note that for K 2 = K 2 = 0, the r. h. s. 
vanishes. Thus in the /C 2 = 0 sector, Gq^ is diagonal. Furthermore, we note that if 
[mi + m\) is odd, the corresponding element of B is zero, signifying that odd and even 
mi values do not mix. 

This, in turn, implies that there are two blocks in B corresponding to every non- 
zero value of /v 2 . One of these, of size [Li/2] (where [m] is the largest integer contained 
in m) corresponds to odd mi and m^ while the other of size [(Li + l)/2] corresponds 
to even values of mi and m^. If we denote these two blocks by B°^[k 2 ) and B [ 1 ^ 2 ) 
(for odd-odd' and “even-even”) respectively, we see that from Eqs.(3.32) and (3.35) 
that their matrix elements, in terms of normalized basis vectors are 


{kiK2\B°°\k,'ik' 2) _ {mi\B°°[K,2)\mfi) 


— (1 — cos K 2 ) cos 
Li 


miTT 


cos 


m'^TT 


(3.36) 


and 

8 1 fO 0'7\ 

{m[\B^%K2)\mi) = j^{l-cosK.2)[l--Smi,o)cos-^cos-^ (. ) 

We point out that since (mi = 0|5^^(«2)|mi) = i(mi|J5“(/C2)|mi = 0), B'^(«2) is not 
symmetric, while B°°[k2) clearly is. 
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It is worth mentioning that the boundary matrix for the finite square lattice com- 
mutes with the corresponding periodic matrix and hence the two can be simultaneously 
diagonalized. Even this turns out to be unnecessary, for that case, because we can make 
B a null operator, satisfying B\k) = 0 \fk, by choosing {x\k) as in Eq.(3.10) and letting 
the origin lie on the dual lattice. 

Returning to Gq^ = P — B, we note that this differs from —B only in its diagonal 
elements. It follows that Gq^ is diagonal in the K 2 = 0 sector and each of the blocks 
in it corresponding to a given k ,2 = n 2 Tr / L 2 ^ 0 can be written as 


^mi ,m2 

-pACm, C'mi+2 + 2,7712 + 2 ^* 071 +2 

P-^C>rrii ^rni -\-q 


-pAC, 


mi+q^mi 




(3.38) 


where is given by substituting values of «i and «2 in Eq.(3.33), p = 1 {00 case), 

p = 1/2 (ee case), A = 8/Li(l - C 0 S/C 2 ) and Cm, = cos(mi7r)/(2Li). The starting 
value of mi is 1 (00 case) or 0 (ee case) while mi+q = Li — 1 {00 case) or Li (ee case). 
Having explicitly determined the forms of P and B in the k- basis we can write 


G„-> = P - (B-” e B”) = © (GJ*)“ (3.39) 

Each of these blocks can be individually diagonalized to arrive at the eigenvalues and 
eigenvectors of Gq^. The eigenbasis of can, unfortunately, not be accessed from 
the position basis through any known fast transform, and we shall, therefore, continue 
working in the k- basis in preference to carry out these diagonalisations. This is 
because much as Go{k,k') is not diagonal, each block in it can be reduced by simple 
manipulations of the rows to an effectively tridiagonal matrix, which can be solved for 
{k\e), in terms of {k\d), in essentially linear order. Since (fe|d) = can 

be obtained from {x\d) through a fast transform and {x\e) = '£k{^\k){k\e) from {k\e) 
likewise, the entire procedure becomes effectively linear. 
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Since is diagonal for /C 2 = 0 one can directly solve for d(nj,n 2 ) instead of 
using the the tridiagonal procedures. Furthermore, in this sector, the restriction of 
d(Ki = 0, K 2 = 0) = 0, removes the zero mode and implies the usage of the gauge- 
jdxing condition 52,- 0, =0. 

To see this in explicit detail we recall that in solving 

G[u] = [u] (3.40) 

where G is of the form given in Eq.(3.38) for [u] in terms of [n], scalar multiplication 
and addition are allowed operations on rows but not on columns[l28]. 

Now, from Eq.(3.38), we see that by multiplying the zeroth row by Cmi+2/(pC'mi) 
and subtracting it from the second, we reduce G to 

'^mi ,m2^mi +2/{]P^mi ) '^mi+2,m2 ^ 

pACmi Gmi-t-2 •^mi+2,m2 ACffii+2GTni+2 

—pAGmiGmi+q ^mi+g,m2 ~ ACmi+qGmi+q _ 

( 3 . 41 ) 

By further multiplying the row by C'mi+2r+2/C'mi+2r and subtracting it from the 
(r + 1)*^ row r = 1, • • • 5 — 1, we turn G into 

-^mi ,m2 ^mi +2/(P^mi ) ■^mi+2,m2 ^ 

0 Ami+2,m2^mi+4/C'mi+2 "-^1111+4, m2 ^ 

-pAC„^C„,+, - AC„.+,C™.+, _ 

( 3 . 42 ) 

i. e. into a matrix which is non-zero only along its diagonal, supen^iiagonal and 
last row. A three-band matrix of this form has (i) the same number of elements as 
tridiagonal matrix and (u) is amenable to a fast 0(N) solution which is similar to that 
used in tridiagonal systems. We remark that tor a given configuration the numbers 
C„.+z.«/(pC„,„.) are constants and need be determined once at the outset. 





We conclude this section by verifying that for a lattice which is rectangularly peri- 
odic in the y-direction, the direct and the inverse transforms are both separable with 
respect to the arguments concerned. More explicitly, 


m=j2{k\x){x\d) 


(3.43) 


= ^d{xi, X 2 )cos ^ 

X ^ 

Tn\KX\ If \ 

COS — 2^ d(xi , X 2 ) 


mx'KXx m2'K{2x2 + 2 : 1 ) 
cos ~ 
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X2 


L 2 

■ ' 2m27rx2 m2'KXi . 2m2'xx2 . m27rxi 

cos — = cos — r sin — sin — = 


where X 2 = 0, 1, ■ • - ^2 — 1 and Xi = 1/2, 3/2, - Li — 1/2. The summation over X 2 can 

be performed by noting that J^xi ^ 2 ) cos and ^ 2 ) sin — are 

the real and imaginary parts of the fourier transform of d(xi, X 2 ). This (fast) fourier 
transform is followed by the multiplication of the real (imaginary) parts by cos 
(— sin ) and the addition of the resulting products, to get the function multiplying 
cos^supr in Eq.(3.43). Finally, the sum over Xi can be performed by using the fast 
cosine transform. The complexity of the whole procedure is then essentially 0{N\nN). 
The backward transform i. e. 


k 


d'(K;i,K2)cos 

mi ,m2 


miTTXi 

~Lr 


cos 


m27r(2x2 4- Xi) 

T 2 


can likewise be performed with a complexity 0{N\a.N). By contrast, for a lattice 
with hexagonal periodicity in the i/-direction, W and W are not separable in Xi and 
X 2 . They can nevertheless be performed with a complexity 0{N In N) using the HDFT 
procedures outlined in[133]. 


3.5 Summsiry and Discussion 

To summarize, we have extended to triangular arrays, the fast algorithms, which were 
so fax available only for square geometries. More specifically, we have developed al- 
gorithms for periodic triangular lattices and those with one set of parallel edges free. 
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In treating the former, we have kept track of the fact that a triangular lattice admits 
rectangular as well as hexagonal periodicity and have shown that the rectangularly pe- 
riodic case is a relatively straight-forward extension of the square case. On the other 
hand, the hexagonally periodic situation is quite distinct in that the eigenfunctions 
of the corresponding discrete laplacian are not separable in Xi and X 2 . Their eigen- 
basis can nevertheless, be accessed from the position basis in 0(N In JV) steps by the 
Hexagonally Discrete Fourier Transforms. This circumstance has allowed us to take 
the algorithm through. The finite case is rather more difficult for the triangular lattice 
than it is for the square. This is because the creation of a boundary in the former, 
requires the removal of two bonds placed obliquely with respect to each other, rather 
than one. This case therefore requires a new approach and this has been developed in 
the paper. In particular we have broken up the inverse laplacian Gq^ into a periodic 
(P) and a boundary (B) matrix and have shown that for a specific choice of the basis 
both B and P — B = Gq^ take on forms which are amenable to diagonalisation in 
linear order. This has allowed us to present for the finite case an algorithm which is 
asymptotically linear to within constants and factors of In N. I would be worthy to 
investigate the case of the fully finite array by removing the y periodicity too. 


We should also mention that we have not fully exploited the symmetries of the 
triangular lattice in this paper. In particular, we still need to investigate the circum- 
stances under which the wavefunctions related by symmetry transformations to the 
one we have actually used (see Eq.(3.23)) are more convenient. 


It is also interesting to note that the approach adopted in this paper has much in 
common with the one used in Sec. 2.4 to treat missing bond defects in square arrays. 
There, too, was decomposed into two parts: one (Gq^) for the perfect finite lattice 
and the other (h) for the defect; G~^ = — /i, with Gq ^ = P—B and [P, B] = 0. For a 

linear defect with n(< Ni) broken bonds, {Go\h] ^ 0. The form of h, however, makes 
the problem amenable to a perturbative approach rather than to the tridiagonalisation 

described in this paper. 
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Finally we observe that for a hexagonally periodic lattice, equivalent points in 
different tiles (see Ref. [133]) can be exchanged to create a new fundamental period. 
These periods therefore come in several shapes of which Fig. 3.1 is but one example. 
The corresponding eigenvalues and eigenfunctions, however, remain unaltered and the 
0{N\nN} algorithm continues to be equally applicable to each of them. These shapes 
can serve as starting points for creating several different finite lattices. The various 
types of boundaries which result are especially significant in view of the interesting 
question; For which class of boundaries does the matrix B commute with P? Indeed 
if there are cases for which [P, P] = 0, they will automatically lead to 0(N In N) 
algorithms. 



Chapter 4 


Defects 


4.1 Introduction 

Recently the study of the dynamical properties of JJAs, especially in the presence of 
controlled [76] and random disorder[74, 109], has attracted lots of attention. In accor- 
dance with Forrester eta7[65], we note that disorder can be introduced into JJAs in 
three ways: (i) Bond disorder in which the coupling strength of individual junctions 
is varied, (ii) Positional disorder in which the vector potential entering via an exter- 
nal magnetic field takes on random values, and (iii) Bond dilution in which certain 
junctions are eliminated altogether. Bond disorder has been studied by Berge[66] to 
distinguish between the Ising and the Kosterlitz-Thouless (KT) transition tempera- 
tures in fully frustrated square lattices and by Chung eta/. [68] and Li et a/. [69] who find 
that disordered arrays have a higher resistance than perfect ones. Positional disorder 
has, similarly, been investigated in various contexts by a number of authors. Reentrant 
glassy behaviour in the temperature versus disorder-parameter plot for positionally- 
disordered arrays has been predicted by Granato eta/. [73] and investigated experimen- 
tally by Benz eta/. [72]. It has further been found that such arrays lead to the formation 
of a novel vortex pattern far from equilibrium [74] as also to a plastic— like flow of vor- 
tices[110]. Finally, experiments on the critical behaviour of randomly bond-diluted 
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arrays[75] have revealed that disorder of this type does not destroy the scale invariance 

I 

of the KT transition. Furthermore, Leath etal.[ll] have discovered that the transition 
to resistive behaviour in arrays with missing junctions occurs only when the external 
current reaches a value for which a vortex path spanning the entire sample becomes 
possible. This is in contrast to breakdown phenomenon in fuse networks and other 
linear systems, wherein breakdown emanates from the most critical defect and spreads 
outwards[131]. They have also found that current-flow past a row of missing bonds 
in an otherwise perfect array gives rise to voltages on becoming resistive. The be- 
haviour of voltages changes from periodic to quasi-periodic and eventually chaotic, as 
the vortex “street”, in the corridor containing the defect, spreads to adjacent columns. 

In this chapter, we investigate additional effects of bond-disorder and dilution in 
current driven JJAs. More specifically, we study, through computer simulations, the 
behaviour of vortices in the superconductive regime, i.e. their nucleation, movement 
and pinning. We do this for arrays much larger than those studied previously. We 
find that the time-independent states of uniformly-driven arrays with linear defects, 
created through bond disorder or dilution, fall naturally into various sectors, separated 
by kink discontinuities in energy. Each sector is characterised by a fixed number of 
pinned vortices. These pinned vortices produce a hysteresis in the system, which we 
explore in some detail. Several complementary insights into each of these phenomena 
are obtained by keeping track of the changes in phase variables caused by modifications 
in current drive. Studying these changes carefully for the periodic regime, in addition 
to the superconductive one, accords us a deeper understanding of how the breakdown 
of superconductive flow actually occurs. We find, for instance, that the periodic flow 
of vortices is sustained by the absorption of energy present in spin-waves. Finally, 
we examine current flow past extended defects. We find that in some cases increased 
bond dilution can enhance the of the system, provided the current flow around the 
extended defect is more streamlined. This effect is analogous to hydrodynamic flows. 

The chapter is organised as follows. In Sec. 4.2, we study bond disorderd/diluted 
linear defects, the transitions to various vortex— sectors, and the periodic regime present 
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in such arrays. In Sec. 4.3, current flow past defects of various shapes is studies to 
explore streamlining properties in analogy with hydrodynamics. The summary and 
conclusions are presented in Sec. 4.4. 


4.2 Linear Defects 

Using the fast algorithms (see Chap. 2), we have numerically studied ^ large arrays (as 
big as 128 x 256 in size) with linear defects of varying sizes. The case of bond disorder 
can easily be incorporated into such a formalism as it affects only the divergence term 
in Eq.(2.4) since iosin^jt giosinOjk. We call such bonds as “graded” bonds. A 
typical geometry is shown in Fig. 4.1. A linear defect, consisting of n broken/graded 
bonds all parallel to the x-axis, is placed symmetrically along the central column (CC). 
The current drive, is uniform along the left edge. We have taken advantage of the 
symmetry, 0{x,y) = 0{x,-y), about the y axis to speed up the simulations 2. We 
thus calculate only the y > 0 section of the array. A positive vortex located at (x, y) 
has a negative image vortex located at (x, — y) due to this symmetry. Henceforth, all 
array and defect sizes refer to the size of the full array and the defect even though 
the calculations are done only on half the array. The number of vortices and their 
description, however, pertains only to y > 0 section of the array. All currents are 

scaled in units of io. 


4.2.1 Graded Defects 

Here we explicitly show that vortices can be pinned inside a linear defect of graded 
bonds. This pinning in turn increases the ig of the system signifying a greater stability. 
We first study the steady-state regime of an array with a linear defect of uniformly 
graded bonds with grading y. The variation of the critical current for the breakdown 


iThe calculations were done on a fiTwlett Packard 9000 series modeipS worksta«on. 

2The r0(x \) = «(-x,y) also exists and could be exploited. However, we have found that 

ine symmetry cirmTnptrieq about the x- and y- axes are broken. 

if the array is driven into the chaotic region, v nnt driven chaotic 

W. monito, both a»d »(-.,!,) to ascertam thot tho art.. It oot dnveo ohaohc. 
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Figure 4.1: Schematic figure of y > 0 portion of a 4 x 16 array with a n = 4 defect 
(shown by a shaded region). Black dots and thick lines joining them symbolize su- 
perconducting islands and Josephson junctions, respectively. Narrow lines inside the 
defect denote graded or broken bonds. Arrows marked i symbolize a current i^xt in the 
direction of the arrow. 

of superconductive flow in the array, with g is shown in Fig. 4.2 for a 64 x 64 array 
with 20 graded bonds forming a linear defect. The most striking feature of the curve is 
minimum at y = Intuitively, one would have expected ic to decrease monotonically 
with g since a smaller value of g implies a reduced capacity to carry super-current. 
However, this is clearly not the case. To understand the origin of this minimum, it is 
useful to examine the behaviour of the phases 0 at a number of relevant sites: 

For very small values of iextt the x-bond current in the CC decreases monotonically 
from a maximum at the center of the defect to ~ i^xt at large values of y. In other 
words, the ^’s along the left edge of the central corridor (LECC) are arranged (see 
Fig. 4.3(a) in a monotonically decreasing sequence (while those along the right edge 
are just their negatives and the phase differences across x-bonds are twice the LECC 
phase values). As i^xt increases, so does each of these phase angles. Eventually, for 
W = ii, the phase, di at the site (-1/2, 1/2) (see Fig. 4.1) reaches a value of 7 r /4 and 
the current sin(20i) in the corresponding x bond becomes unity, i.e., it turns critical. 
As iext is raised beyond ii, the phases simply go on increasing. As a consequence, 2di 
goes into the second quadrant and the current in the junction begins to fall. Thus, 
the criticality of a single junction does not mark the onset of dissipation, as also noted 
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Figure 4.2: A plot of ic and i„ vs. g showing the minimum in ic for 64 x 64 array with 
a linear defect consisting of 20 graded bonds. 
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by Leath and Xia[77]. However, in contrast to their observations, this criticality is not 
accompanied by the formation of a vortex. A further increase in i^xt causes the next 
few bonds in the corridor away from the center to successively carry the critical current. 
This phenomenon persists until i^xt reaches a value di becomes equal to 7r/2 and the 
current in the junction drops to zero. As soon as iext exceeds the current becomes 
negative and a vortex appears since the next x-bond still carries a positive current and 
there is a current circulation in the central plaquette located at (0,1). The behaviour 
of this vortex now depends on whether g > Qm ot g < g^- For 0 < ^ < 5r,„, the vortex 
remains pinned inside the defect (Fig. 4.3(b) with the pinning center shifting to larger 
values of y as g increases. Eventually, aA. g — g^, the pinning occurs just outside the 
defect at its tip. For g > gm, there is no pinning at all. We notice, however, that the 
vortex is never pinned at its point of nucleation even as g ff. 

Finally, for g < gm but < i^xt < ic-, the pinned vortex moves outwards, only to 
get pinned again, always inside, or at most at the edge of the defect. At ie^t = 
the Lorentz force on a vortex at any point along the CC exceeds the pinning force. 
Consequently, the vortex moves all the way out to the edge of the array and begins to 
cause dissipation. 

The i’u vs. g curve (see Fig. 4.2), on the other hand, does not show any minimum. 
For g > gm the and the curves are coincident but separate out for g < gm- We 
thus infer that grading causes the pinning of vortices which in turn raises the ic of 
such systems. This is quite in conformity with the case of continuum superconductors 
where vortex pinning centers are artificially created (e.g. by ion bombardment [147]) 
to avoid dissipation. What is surprising, however, is that gm ~ 0.2 is independent of 
the lattice and defect size provided the latter is small compared to Ny. In the case of 
ra ~ iVy the vortices do not get pinned for any value of g due to the presence of edge 
effects in the array. Consequently, for such an array increases monotonically with g 
without showing any minimum. 

The profile of currents in the x-bonds in the CC, for the case described above, 
suggest that it is possible to pin several vortices inside the defect for g < gm provided 
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(a) (b) (c) 


Figure 4.3: Typical configuration of phases in the central corridor of a 32 x 128 array 
with a n = 50, fir = 0.1 graded defect. Only a portion of the corridor is shown. The 
shaded region is the defect. The vortices are marked with a + sign. For (a) i^xt < iv 
while for (b) and (c) < i^xt < 4- The number of pinned vortices in (a),(b) and (c) 

are 0, 1 and 2 respectively. 



70 


the latter is sufficiently long. In accordance with this expectation, we have found that 
for a long narrow array 8 x 128 with n = 60 (g = 0.05) and for iext — 0.064 we get as 
many as 3 vortices pinned inside the defect (as in Fig. 4.3(c) which shows the presence 
of 2 vortices). After the formation and the movement of the first vortex, the magnitude 
of 6i at the center of the array once again passes through the same cycle of events to 
produce a second vortex which forces the first one outward (the second vortex is formed 
when 6i = 37r/2). From a plot of the x bond-current in the CC versus y, one can easily 
locate the pinned vortices since these are clearly situated wherever a change in sign of 
the bond-current occurs. 

4.2.2 Broken Bonds 

We next investigate the behaviour of arrays in the presence of a linear defect produced 
by bond removal. It had been conjectured[77] that for infinite lattices, the pinning of 
vortices nucleated by a defect of this type, would occur for defect sizes 200. With 
our improved algorithm, we are in a position to check this out. 

To begin with, we observe that the sequence of events described in Sec.(4.2.1) is 
followed for this system as well, with some minor differences: The vortices are now 
formed at the tip, rather than the center, of the defect and their subsequent behaviour 
depends on the value of as compared to ip, the depinning current required to overcome 
the potential barrier[16] for vortex-motion from one plaquette to another. If < ip, 
the vortices are pinned by the lattice away from the defect. The transition to the 
periodic (resistive) state takes place at i^xt — '^c ^ iv For i„ > ip, the vortices move 
right across the central corridor, out of the array, and the transition to the periodic 
state takes place without any vortex-pinning. In this case, i^ coincides with i„. 

We next describe how ip, i„, and i^ vary with Nx, Ny, and n. This dependence 
leads to an understanding not only of the Nx,Ny-^oo limit, but also of a number of 
finite-size effects which are interesting in their own right. 

We recall that Lobb et a/. [16] calculated the potential barrier, Eb for a vortex to 
move from the center of one plaquette to the adjacent one for square arrays assuming 
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a sinusoidal form for the pinning potential. Their fonnalism[l48] can likewise be used 
to calculate ip = Eb/(2Ej) for various values of N^,Ny, where Ej is the Josephson 
coupling energy. In carrying this out for vortex motion in the j/-direction, we find 
that for fixed N,=32 and = 4, 8, 16, 32, i, = 0.147, 0.112, 0.102, 0.099, respectively. 
Similarly for AT,, = 32 and = 8, 16, ip = 0.072, 0.093, respectively. thus see that 
in rectangular arrays, the finite size effects in the x— and y— directions are different, 
resulting in a highly anisotropic ip, which can, for < Ny, be larger than the infinite 
array limit of i^ ~ 0.1[148]. This implies that in finite arrays, it should be possible to 
observe vortex pinning for a defect of size smaller than ~ 200 [77]. This is confirmed 
by our simulations, e.g., of a narrow 4 x 256 array with n = 22. Here we observe as 
many as 10 pinned vortices lined up in the central corridor for iext = 0.17396, which 
incidently is much larger than 

To isolate the effect of finite array size on i,.? it is expedient to define Ai^ — — ic^ 

where if is the critical current in an infinite array. A log-log plot of Aid if vs. Nxln 
for a fixed Ay/n = 32 is shown in Fig. 4.4. In obtaining this curve, we have taken 
if = ic{Nx = 128) for all values of n considered. To convince ourselves that an array 
with Nx — 128 is infinite for all practical purposes, we have checked that i^ varies by 
less than 0.3% and 0.002% as Nx is changed from 64 to 128 for n = 8 and n = 2, 
respectively. From Fig. 4.4, we see that 4 scales in Ar/n and that, furthermore, for 
n ~ 2Aj,, Aid if is ^ 30%. This can be understood as follows. The current injected 
at the left edge of the array must flow around the defect. The smaller the value of Aa,, 
for a given defect size, the larger the curvature of the current flow lines around the 
defect and hence the larger the vorticity dixjdy. Since ic is related to nucleation of 
vortices, and this occurs more easily for an enhanced vorticity, 4 diminishes, if Nx is 
reduced keeping n fixed. 

A similar reduction in ic also occurs with decreasing Ay for a given Nx/n because 
current redistribution away from the defect gets increasingly constrained. In Fig. 4.4, 
the error in ic due to the finiteness of Ay is ~ 0.01%. 
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Figure 4.4; A log-log plot of (Aic/i“) vs. (Nx/n) for various n but constant Ny/n = 32 
showing the scaling behavior. 


The effects of a sharp increase in vorticity near the tip of the defect are dramatically 
exhibited by the formation and flow of vortices along rather than perpendicular to the 
direction of i^xt- This is observed for small values of A^/n for which the current flowing 
near the tip has a large dix/dy component. This in turn produces large values of both 
iy and diyjdx. As a result, pairs of vortices ( both of the same sign due to symmetry 
about the y— axis) nucleate at the sides of the tip for i > and flow in opposite 
directions along the x-axis. This phenomenon can be observed, for example, in an 
8 X 128 array with n = 48 and i^xt = 0.15. 

Having shown that in long narrow arrays with moderate-sized defects, ip goes up 
and outstrips while in general comes down, we turn to a detailed study of the 
vortex-pinning regime > ip > i^xt > iv)- For concreteness, we specialise all our 
observations to a 16 x 256 array with 100 missing bonds. We find that the vortex 
which appears at the tip of the defect for text = iv is always unstable there and moves 
up the corridor as soon as it forms. However, it gets pinned 15 plaquettes away from 
the defect and all voltages then quickly drop to zero. The phases along the LECC are 
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now configured as follows: Those preceding the vortex are all in the second and third 
quadrants while those following it are in a decreasing sequence lying entirely in the 
first quadrant. A continuous enhancement in sees each of these angles increasing 
continuously. Whenever the phase difference across the upper edge of the plaquette 
containing the vortex becomes tt, the current through it changes from being parallel 
to iext to being antiparallel and the vortex moves up by one plaquette. Likewise, when 
2etip reaches the value {2k + l)7r, fc = 1, 2, . . ., the vortex is nucleated at the tip of 
the defect. This vortex immediately moves up the corridor pushing the other vortices 
upward ahead of itself. All the vortices which are thus set into motion eventually get 
pinned and the transient time-dependence in voltage disappears. The process continues 
until (for the above-mentioned values of Nx, Ny and n) a maximum of five vortices are 
pinned by the lattice. Thereafter, for iext > io the vortex street begins to run with a 
definite periodicity. It is noteworthy that in the periodic domain, the minimum number 
of vortices simultaneously present in the central corridor either equals or exceeds by 
one the maximum number pinned in the steady-state regime. 

The vortex train formed in the CC is found to be unstable in long narrow arrays at 
any finite temperature (the rounding off errors in numerical simulations are equivalent 
to a small but finite temperature). For long narrow arrays (e.g. 4 x 32), the vortices 
in the CC are close together, typically 3-4 lattice spacings apart. This is because the 
image charges across the edges in the narrow dimension screen the vortex charges in 
the CC. As the position of a given vortex inside its plaquette fluctuates away from 
the center, a; = 0, it experiences a destabilising force, due to y bond currents in 
this plaquette, produced by the other vortices in the CC. This current is large enough 
(owing to the proximity of these vortices) to overcome the pinning force in the x— 
direction. The images of the vortex across the width of the array are also sufficiently 
close to exert a force which further helps the destabilisation process. The linear chain 
of the vortices in the CC is thus rendered unstable. This instability is also seen in 
wider arrays with long linear defects, which create appreciable y-bond currents in the 
plaquettes along the CC close to the tip of the defect. 
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Figure 4.5: A plot of critical currents (•) and ic {x) vs. n for a 16 x 256 array. 
For n ~ 74,2c = 2 „. Also shown are if* values (o) using corrections for finite Nx from 
Fig. 3. The inset shows a magnified region where and ic bifurcate. 
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Finally in this section, we consider the asymptotic dependence of on the defect 
size n, particularly in the presence of pinned vortices. Since vortices redistribute the 
current, increasing its value away from the defect, they could conceivably make the 
defect effectively larger[77]. On the contrary we find that, as far as is concerned, the 
effective size of the defect is reduced by pinned vortices. We plot vs. n in Fig. 4.5 
for a 16 X 256 array. For n ^ 74, the vortices cannot be pinned as Zc(= iv) is larger 
than ip, while for n ^ 74, they can. Hence, for n > 74, we plot both (full circles) 
and ic (crosses). It is seen that while the values of follow the curve smoothly across 
the value n = 74, those of ic bifurcate from and lie above this curve. This increase 
in ic relative to iy shows that the effective size of the defect is reduced by the pinned 
vortices. The extrapolated values of if’ using the results of Fig. 4.4 are also shown. 
Each of these values is larger than i^(c^ 0.1). It follows that defects of size much larger 
than 74 would be needed to observe pinning in infinite arrays: The estimate of Leath 
and Xia[77] is in the right ball-park. 

4.2.3 Energy Of the System and Hysteresis 

The energies of the static configurations forming the steady-state regime of a bond- 
diluted array acquire, in the presence of pinned vortices, a range of interesting features 
which we now investigate. As mentioned in Sec. (4.2.2), a vortex which forms at i^^t == '^v 
never stabilises at the point of nucleation (the tip of the defect for the case of broken 
bonds) and moves a certain distance along the CC before getting pinning. Thus the 
variation of the phases in the CC (and elsewhere) is discontinuous between and i+. 
It follows that the energy of the system defined as[62] 

= V [1 - cos %] - ie.t E 

<ti> 

registers a downward discontinuity at the formation and eventual pinning of each new 
vortex. The energy versus w graph (see Fig. 4.6 ) can thus be conveniently divided 
into different sectors with the sector having p pinned vortices. The energy decreases 
continuously within a given vortex sector and registers a kink discontinuity only where 
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Figure 4.6: The energy of 16 x 256 array with n = 100 the system as a function of 
iext < ic- Discontinuities appear at each transition from one vortex sector to the next. 
The inset shows a magnified transition region. 


one enters a higher vortex sector. Moreover, the difference between successive values of 
text, at which these kinks occur, decreases continuously, while the energy discontinuity 
itself becomes larger and larger. 

These observations suggest that for an infinite array (or one which is infinite in 
the y-direction, at least), there would be an infinite number of vortex sectors. The 
energy of the system for any of these sectors would be infinite, but the energy differences 
would be finite and well-defined. However, the sequence En of energy-gap values would 
diverge for n — »• oo. On the other hand, the sequence of ie^f-values marking transitions 
firom one vortex sector to another would rapidly converge to the depinning current, 
i'^ 1 for a single vortex in a perfect array. The spacing between vortices pinned by 
the lattice would increase with distance from the defect and would become infinite for 
those furthest out. The vortex— vortex interaction, which inevitably comes into play 
for finite separation between vortices, would thus become negligible at large distances 
from the defect[126]. It is because of this that the infinite vortex street would be set 
into motion only by iff. 
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Figure 4.7; The hysteresis loop for a 8 x 256 array with n = 50 for iext < ^c- 
figure the y scale represents energy/bond. 


Finally, once vortices have formed and stabilised at various pinning sites in the 
presence of a certain iexu they will remain as configured even when this current is 
reduced to zero. There are consequently at least two configurations of the array at 
iext = 0 with different values of energy. The actual number is, of course, much larger. 
In fact there are as many states as there are allowed vortex configurations in the time- 
independent regime of the given array. The energy is, therefore, multiple-valued in 
iext and the curve is, in this sense, hysteretic. A typical graph, for an array 

of size 8 X 256 with 50 broken bonds, is shown in Fig. 4.7. We note that this graph 
is symmetric about the E-axis. In plotting it, we first raise iext from 0 to t,nax, a 
current for which there are 5 pinned vortices in the array, and obtain section OA of the 
hysteresis curve. We next decrement the current from through 0 to -t^^x (curve 
APQBC) and finally raise it back to w (curve CPRDA). The hysteresis loop so 
obtained is found to be independent of the size of the current-step used in tracing it. 
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The various features of the closed curve displayed in Fig. 4.7 are understood as 
follows. A set of 5 vortices forms and stabilises over part OA, as described above. There 
is no change in the positions of these vortices as the array goes from A to Q. However, 
the system’s energy increases steeply (so much so that at text — 0) corresponding to 
point P in Fig. 4.7, the energy is positive). Both circumstances persist until we come 
close to <5? wherein one of the vortices gets annihilated by the anti-vortex nucleated 
at the tip of the defect decreasing its energy. Slightly short of this point, anti- vortices 
steadily begin nucleating at the tip of the defect. Once this happens, the vortices 
originally present in the CC rapidly disappear through vortex-anti- vortex annihilation 
and their positions are taken shortly thereafter by anti-vortices. The energy of the 
array suffers a sharp drop in the process. For i^xt — —imaxi this energy is negative and 
equal to its value at imax- Since the figure is symmetric about iext = 0, a similar process 
occurs during the branch CPRDA with the vortices now annihilating the anti- vortices 
present in the array. Interestingly enough, at our level of resolution in iext > we did not 
encounter, at any point of the hysteresis loop, a regime in which there are no vortices 
present in the array. Finally, it should be pointed out that hysteresis loops of this kind 
are not uncommon among superconducting systems [149] 

4.2.4 The Periodic Regime 

In this section we examine the periodic regime of arrays with linear defects and study 
the mechanism of propagation of vortices in the CC. Once the system enters the time- 
dependent state (for iext > Q, vortices nucleate and peel off periodically from the tip 
of the defect, traverse the entire length of the corridor and pass out of the array [77]. 
All the phases in the array rotate at the same average angular velocity, u, where the 
averaging is done over a time period. The angular velocity is not constant since that 
would imply that the super-current in each link is independent of time. 

A snapshot of the phases in an array of size 32 x 128 is shown in Fig. 4.8. All 
the phases to the left of the CC rotate anticlockwise while those on the right rotate 
clockwise. The phases along the columns of the array form “spin waves” moving in the 
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Figure 4.8: A snapshot of an 32 x 128 with n = 9 array in the periodic reginae. The 
arrows represent the phases at the superconducting sites. The defect is the lightly 
shaded region. The dark shaded plaquettes are occupied by vortices. The wavefronts 
corresponding to the spin waves are seen as dark streaks across the picture and form a 
wake in front of the vortices. Note that the wake is asymmetric about the y-axis. 



positive y— direction, while the phases along the rows constitute spin waves converging 
onto the CC. The resultant spin waves, therefore, seemingly form a wake, but actually 
move into the vortex and propel it forward. In Fig. 4.8, the wavefronts corresponding to 
the spin waves can be seen easily as dark streaks running across the picture. However, 
these wavefronts are ahead of the vortex rather than behind it. The situation is akin 
to a time reversed motion of a boat moving in water. Whereas the boat provides 
the energy to propel itself forward while transferring some of it to the waves in the 
wake, here the vortex absorbs energy to move forward from the spin waves present in 
the JJA. The “wake” here is asymmetric along the x-direction as the vortex has the 
special property of changing the phase by tt across the plaquette containing it. 

The nature of the waves discussed here is different from those in arrays with ca- 
pacitive junctions[53, 136]. In such arrays, a vortex moving across a plaquette induces 
voltages on the superconducting islands, which keep oscillating due to stored capac- 
itive energy long after the vortex has moved away. Thus, vortices (having a finite 
mass in capacitive arrays) can interact with each other through these spin waves. In 
overdamped arrays, the vortices just “ride” the spin waves described in this paper. In 
capacitive arrays, we observe that the inertial oscillations are superimposed upon the 
waves that the vortices ride. 

The separation between successive vortices along the CC is A/2, where A is the 
wavelength for spin waves in that direction. The frequency u) is roughly proportional 
to 1/A showing a linear dispersion. The value of A is, however, affected by edges. 
It increases monotonically as one moves towards the boundary. This is because the 
frequency of the periodic regime is constant at lo (for a given while the speed of 

vortices is higher near the edges due to attraction from the image charges across the 
latter. 

In the x-direction, the phases wind with a finite wavelength even in the steady 
state regime. This static configuration corresponds to a wavelength which decreases 
monotonically with an increase in i^^t- As i^xt crosses this “frozen” spin wave starts 
moving. Consequently, at = i,, w = 0, but the wavelength in the x-direction is 
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finite. Sinailarly a finite wavelength along the ^direction would also exist for u; = 0 
in the case of sufficiently large arrays where pinned vortices in the CC are present in 
the steady state regime (and the transition from steady state to periodicity takes place 
through an n-vortex sector with finite n, as described in Sec. 4.2.3 above). 

The angle subtended by the wake at each vortex can be deduced from the ratio of 
the wavelengths in the y— and x— directions. As expected, this angle decreases as the 
vortices move faster under an increase in 


4.3 Extended Defects 

We now turn to more complicated defect-patterns with a view to determining the 
dependence of ig on defect-shape and size and distribution. In doing this, we shall 
find it useful to keep in mind the analogies, brought out by Mehrotra and Shenoy[78, 
79], between current flow in JJAs and hydrodynamic flow through a pipe. These 
authors consider non-uniformly driven JJAs at zero-temperature. They find that the 
linear gradient of the external current drive, dix/dy is a measure of injected vorticity 
which, multiplied by the y-dimension of the array (pipe diameter), and scaled by 
the inverse shunt resistance gives the analog of the Reynold’s number. They further 
note that the chaotic behaviour observed above a certain injected-vorticity threshold 
is, as in turbulent fluid flows, due to a mixing of positive and negative vortices. In 
the spirit of this analogy, we now liken current flows past bond diluted/disordered 
linear defects to fluid flows past obstacles. Such current flows are closely related to 
those explored in Ref. [79] because the defect automatically creates a non— zero current 
gradient, {dQ /(dy). We examine effects of streamlining and find that similarities with 

fluid flow continue to exist. 

We study the bond diluted case by grading defective bonds with a y 0 (we use 
g = 10-’’, in actual practice) rather than by eliminating them. This permits the use of 
faster algorithms and makes no difference to the ic valines we are interested in. Indeed, 
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Figure 4.9: The critical current vs. the separation d between defects. 

only supercurrents flow through the bonds in the tinae independent state, and even 
these turn olf wherever g is set to zero. 

To determine the effect on ic of the distance between defects, we consider the case of 
two highly graded bonds [g = 10“^), placed in the central column with an inter-defect 
separation d. A plot of versus d (Fig. 4.9) shows a minimum at d = 2. As both 
defects try to force current away from themselves, a “squeezing” of current into the 
channel between the two defects occurs resulting in a “breakdown” inside the channel. 
This IS similar to hydrodynamic flow through a narrow passage. For d < 2, the channel 
is too narrow for current to flow through and the two defects are seen as just one large 
defect. This causes an initial decrease in i,'. For d > 2, the current starts flowing 
through the channel. For large values of d, the two defects become independent of each 
other and approaches its value for a single defect. 

In studying the effects of streamlining, we use the following notation to denote defect 
configurations. A sequence of numbers a,m-b,c,m-b,a corresponds to a defect placed 
at the center of the array, with c graded i-bonds in the central corridor, m-columns 
of b graded r-bonds each followed by a single column of a graded r-bonds, on either 




CHAPTER 4. DEFECTf^ 


83 



Figure 4.10: Schematic figure of the y > 0 section an array with a typical defect 
(shaded) labeled 2, 3 — 1, 3, 3 — 1, 2 used to study streamlining. Each side of dotted 
squares denotes a Josephson junction. Current is injected/withdrawn from left/right 
edge. All the junctions inside the shaded region are graded with g = 10“'^ and are 
equivalent to broken bonds. 

side of the CC. A typical defect configuration (2, 3 — 1, 3, 3 — 1, 2) is shown in the inset 
of Fig. 4.10. 

Our results are shown in Fig. 4.11 where the values of ic are plotted versus m for 
two distinct cases. We first consider the 2,m — 0,3,?ti — 0,2 class of defects which 
corresponds to placing an 2 defective bonds in front of 3, along the injected current 
path. For m. = 0, the ic value' is 0.592, which is 0.025 higher than 4^^ the value for a 
single n = 3 defect. As m is increased, ic goes through a maximum at m = rrimax — 2 
and then decreases asymptotically to 4^^. A second class of defects 2, m — 1,3, m — 1,2 
shows a similar trend, except that the asymptotic value now corresponds to the value 
for the oo - 1, 3, CO - 1 defect (The dotted line actually represents measurements on 
al5-l,3,15-l defect). It should be noted that while the second class has more 
defective bonds than the first, the ic values of the former are higher than those of the 
latter, for all m. 

The above results are consistent with hydrodynamic flows and can be understood 
as follows: Any reduction in the curvature and y-gradient of current density in the 
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Figure 4.11: The variation of ic with m used to study streamlining. The circles rep- 
resent ic values for defects 2, m - 0, 3, m - 0, 2 while crosses are for defect pattern 
2, m - 1, 3, , m - 1, 2. The for single n = 3 defect is shown by a bold line while the 
dotted line is for a 15 — 1, 3, 15 — 1 defect. 


x-direction {di^ldy) increases the value of ic, as it corresponds to an improved stream- 
lining of the current flow. For example, the ic value for m = 0 is higher than i^^^ 
Similarly, rectangular defects have higher ic’s when aligned along the current direction 
as opposed to being perpendicular[104]. As m increases to m^a®, ic increases further 
because of even better streamlining. The current flow does not sense the small sepa- 
ration between the n = 2 and n = 3 defects and does not enter the inter-defect region 
significantly. For m > rrimax, the current starts penetrating this region and in going 
around the n = 3 defect produces a larger value of dixidy. This in turn decreases the 
values of ic- For very large values of m, the defects becomes independent and hence 
ic asymptotically approaches 4^^. As the current flow lines cannot bend back into the 
region between the n = 3 and n = 2 defects far enough in the second class of defects, 
the ic values are higher than the corresponding ones in the first class, for all m. 
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4.4 Summary and Conclusions 

In conclusion, we have systematically studied JJAs with defects. Using a fast algorithm 
we have investigated the steady-state regime, the pinning of vortices, the hysteresis 
present in such arrays and the periodic regime for the case of bond diluted/disordered 
defects. For bond disordered arrays with linear defects we have shown that for g < 
vortex pining increases its stability against breakdown. In the case of bond diluted 
linear defects, we have shown that (i) scales with iV„/n (for fixed Ny/n), (ii) vortices 
can be pinned in smaller arrays than thought earlier, (iii) for narrow (small A^/n) 
arrays the vorticity at the tip of the defect can increase to a point where vortices form 
and move along the injected current direction, (iv) the presence of pinned vortices 
increases the super-current carrying capacity of the array, (v) the energy of these 
arrays registers a discontinuous drop at each transition from one vortex sector to a 
higher one, and (v) these pinned vortices cause hysteresis. An investigation of the 
periodic regime of such arrays shows the presence of waves which propel the vortices 
forward. The study of extended defects have shown that an improved streamlining can 
enhance the ic of the defects analogous to hydrodynamic flows. 



Chapter 5 

Ballistic Vortices 


The dynamics of vortices in two-dimensional JJAs has been a topic of much current 
interest both theoretically and experimentally. A large part of this interest has focused 
on the motion of vortices in JJA and on understanding the physical content of an ef- 
fective equation of their motion[53, 54, 106, 136, 150]. An exciting possibility to result 
from these studies is that of vortices moving ballistically in arrays of underdamped 
junctions[54, 151]. The first experiment which claimed to have observed such ballistic 
motion in the classical regime was reported by van der Zant et a/. [3] (hereafter referred 
to as ZFOM). Computer simulations[56-59, 137], however, failed to observe any bal- 
listic motion in square or triangular arrays. In these simulations a single vortex was 
accelerated by an applied external current which was subsequently switched off. A bal- 
listic vortex would have continued moving under its own inertia. Several explanations 
have been put forward to account for the discrepancy 

Eckern et al.[53] employed a path— integral approach to arrive at an action which 
included the quasi-particle tunneling present in capacitive junctions. Separating the 
field variables into two parts (one of which is due to a single vortex configuration) they 
obtain an equation of motion of the vortex coordinate in the classical regime. This 
formulation leads to a finite vortex mass proportional to the inter-junction capacitance 
of the JJA. Later Eckern et al.[136] show that the vortex velocity is strongly affected 
by the pinning potential of the lattice. As a result they predict that the window for 
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ballistic motion in a square array is negligibly small The case of the triangular array 
is presented as a more promising candidate since the potential barrier was ~ 5 times 
lower in such arrays. 

Bobbert[56], using a piecewise linear potential instead of the sinusoidal pinning 
potential, found a non- vanishing viscosity for low damping and an energy transfer by 
the moving vortex to junctions just behind it. He could not find ballistic motion in 
square and triangular arrays and conjectured that the experiments could be measur- 
ing some collective property involving vortex— vortex interactions (VVI) rather than 
ballistic motion. 

Geigenmiiller et aZ. [57] probed the extra losses resulting from the generation of 
charge oscillations by the moving vortex to the spin waves. They found that the friction 
coefficient rj is proportional to (3^'^ and was independent of the resistance R for large 
values of /?c- Using a continuum theory in the presence of a moving vortex (represented 
by a star-like configuration of the phases) they calculated the vortex velocity as a 
function of the external current. They concluded that ballistic motion could exist only 
under very special circumstances : a weak pinning potential (a triangular as opposed to 
a square array), a driving current just above the depinning threshold ip, and extremely 
high values of Pc- 

Yu and Stroud et aZ. [58] studied tl^e I-V characteristics of an underdamped array 
with / near 0 and 1/2. They found that there exists two critical currents Iciif) 3^d 
Ic 2 {f) for JJAs. While the single vortex depinned at the former value, the whole 
vortex lattice depinned at the latter. The motion of the vortex was non-hysteretic 
although single junctions were highly underdamped. For high Pc and high vortex 
velocities row switching occurred instead of ballistic motion. They also confirmed[59] 
and extended the predictions of Geigenmiiller eZ aZ. [57] to triangular arrays that there 
is extra damping due to loss of vortex energy to the plasma oscillations. 

Fazio et aZ.[60] took the effects of the discrete nature of charges into account the- 
oretically and concluded that the spin wave spectrum is altered so as to reduce the 
dissipation. Ballistic motion of vortices thus becomes possible for a wider range of 
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velocities. Recently, Hagenaars et a/. [137] numerically studied the enhanced nonlinear 

i 

viscosity for vortex motion. They found that the viscosity varies as 77 = yl/(l + Bv"^) 
where A and B are constants and v is the velocity of the vortex. They, however, did 
not find any ballistic motion and concluded that further theoretical and experimental 
work was needed to resolve the discrepancy between experiment and theory. 

In this chapter, we explicitly demonstrate that a vortex-street can move through a 
field-free region, even for a square array with overdamped junctions. In other words, 
we show that VVI are large enough to overcome the pinning potential for vortices and 
cause them to move without any help from an external current drive. Our results thus 
explain the experimental observations of ZFOM in terms of VVI. 

5.1 System Configuration 

We recall that the experimental arrangement of ZFOM involves a H shaped array of 
triangular Josephson junctions. The setup can be divided into three regions ; 

(i) The accelerator region where vortices, produced by an applied magnetic field, are 
accelerated by an external current ig^t- 

(ii) These vortices are collimated by a narrow channel and 

(iii) A field-free detector region to detect the collimated vortices at the far end of the 
accelerator. 

Voltage probes at various points detect the movement of the vortex through the array. 

We have simplified the ZFOM geometry to the one shown in Fig. 5.1. A linear 
defect, consisting of n broken bonds parallel to the x-axis, is placed in the central 
column (CC) at the bottom edge of an N, x Ny square array. Most of our work is 
on two sample arrays, A and B, 16 x 64 and 32 x 128 in size with n = 15 and 30, 
respectively. In sample A(B) a d.c. drive is applied to rows 1-25 (1-45) and 
is linearly decreased to zero between rows 25-35 (45-55) to avoid large gradients of 
current in the y-direction which would cause movement of vortices in the x-direction 
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Figure 5.1: A schematic of a 4 x 8 array with a n = 2 defect. Black dots and lines 
joining them symbolize superconducting islands and Josephson junctions, respectively. 
Arrows marked i symbolize a current in the direction of the arrow. 


as described in Chap. 4. We do not require a channel to collimate the vortices as 
these nucleate at the tip of the defect and automatically move exclusively along the 
CC[76, 77], for the range of currents used. The bus bars used by ZFOM to provide 
a parabolic potential for confining vortices near the CC are likewise not needed and 
we consequently use free boundary conditions. Moreover, by placing the defect at 
the bottom edge, vortices of a particular sign are created as in the ZFOM geometry. 
Finally, for sample A{B), the rows > 41(72) are essentially field-free and correspond 
to the detector region. The current flowing across the CC x-bonds falls to O.liext 
rows 41(72) and to O.Oliext at rows 52(95), as measured in the absence of vortices. 

We emphasize that the arrangement used by us is very convenient to study the 
issue at hand since it automatically produces a collimated vortex beam. The geometry 
is also easily realizable experimentally. 
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5.2 Numerical Simulation 

We have numerically studied the above-mentioned configuration using the standard 
classical RCSJ model for 0 < /3c < 100 and temperature 0 <T < O.OlEj (introduced 
via a Langevin noise term), where Ej is the Josephson coupling energy. We scale all 
currents by io, the single junction critical current, and set the lattice constant a = 1. 
Time is measured in units of = h/{2eRio), the inverse characteristic Josephson 
frequency or Ihe inverse plasma frequency. We use a fast cosine 

transform algorithm extended to arrays with missing bonds as discussed in Chap. 2. 

For 0.25 ^ iext ~ 0.4, which is large enough to nucleate vortices at the defect but 
insufficient to cause column switching, we observe that a vortex nucleated at the tip of 
the defect and accelerated by iext moves up the CC, where it gets pinned in the field- 
depleted region. When another vortex is nucleated, the first vortex is pushed forward 
and both vortices move until each one gets pinned in the CC. This process continues 
until in the s-vortex sector (see Chap. 4) (s is between 7 and 8 for array A), the vortex 
farthest from the defect is pushed out of the array. At this stage the vortex-street 
begins to run and vortices flow out of the array periodically now moving all the way 
across the field-free region to the far end of the array for all the above values of /3c. 

To make direct contact with the experiment of ZFOM, we measure in array A the 
voltage between the current driven edges of the array, averaged both along the edges 
from rows 1—39 and in time. This corresponds to the voltage in the accelerator region 
of ZFOM. Similarly, the average voltages Veh and Vj, from rows 40-41 and rows 42-64, 
correspond to the channel and the detector regions, respectively. We also monitor the 
voltage Ve across the bond in the CC at the far end where the vortices leave the array. 
In Fig. 5.2, we plot each of these voltages as functions of iext at two values of T with 
Pc — 10. For T = 0.003, all voltages are equal (Vch not shown for clarity) showing 
that all the vortices nucleated at the tip of the defect manage to reach the edge of 
the detector region and do not diffuse away from the CC. At a higher temperature, 
T — 0.01, the voltages fall as the vortices move through the accelerator, channel and 
detector regions showing that not all vortices reach the detector region. Furthermore, 
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Figure 5.2: A plot of average voltages in various regions of the array versus current igxt 
applied in the accelerator region. 

the voltage Vf. is smaller than Vj showing that the vortex beam spreads and not all 
the vortices which hit the far end of the array do so right in the center. This behavior 
is just like the one seen by ZFOM. However, the typical value of the temperature Td 
at which vortices start diffusing is found to be ~ O.OlEj as compared to ~ O.lEj by 
ZFOM. 

There are three reasons for this discrepancy: The pinning potential in a square 
array (our case) is ~ 5 times larger than for a triangular array (ZFOM). With forces 
between vortices varying inversely as the separation between them, the linear density 
of vortices for a square array is ~ 5 times larger than that for a triangular array. At a 
given temperature with a certain density of thermally created vortices the annihilation 
of vortices in the vortex-street is then ~ 5 times larger in a square array (cf. law of mass 
action). The loss of a single vortex from the vortex— street disturbs the WI mediated 
motion of the whole vortex-street. Secondly, we have no bus bars and any vortex in 
the vortex-street which gets out of the CC easily diffuses towards the current-driven 
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edges. Thirdly, ZFOM have reduced the io of junctions in the channel and detector 
region to facilitate the entry of vortices in that region. We find in our simulations that 
having bus bars and changing the io of selected junctions increases the ratio at 

a given T. All these effects are of the correct sign and order of magnitude to roughly 
give agreement with ZFOM for the values of Td- 


5.3 Vortex Vortex Interactions 

To trace the origin of the propagation of vortices in the field-free region we first study 
quantitatively the interaction between vortices in finite-size arrays. We take boundaries 
into account by introducing image charges across each edge of the array. We check all 
our calculations against simulations in arrays of sizes 16 x 128 and 32 x 64. A uniform 
current drive of 0.9 initially applied to all rows is turned off after a desired number of 
vortices have been injected into the array from a defect with n = 5. 

We first consider the case of a single injected vortex. This gets pinned in the CC, 
all voltages drop to zero and only supercurrents due to the vortex flow after i^xt is 
turned off. This pinned vortex can be moved by applying a local current id to each of 
the two left corners of the plaquette containing the vortex and withdrawing it from the 
two right corners. With the vortex positioned close to the center of the array, the effect 
of its image charges across the edges is negligible. We now vary id and measure the 
threshold current it required to depin and move the vortex to an adjacent plaquette 
in the CC. We find that it = 0.167, independent of the direction of depinning. The 
effective current exerting a Lorentz force on the vortex is the average current flowing in 
the top and bottom bonds of the plaquette containing the vortex. This can be fotind 
by measuring the fraction a of the current id which flows through these bonds in the 
absence of the vortex. We measure a = 0.6338. Hence, the depinning current ip for a 

single vortex is given by ip = ait = 0.106, in agreement with the predicted theoretical 
value[61]. 
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Figure 5.3: The images of the charge at 0 in the x direction and the test charge P 

We consider next the VVI between two vortices located at y, and yj in the central 
column (point O and P in Fig. 5.3 respectively). As topological charges in 2D they 
interact via a logarithmic potential, — 27rEj-logy,j where = \yi — %|. The repulsive 
force, Fij(yi,yj), between the vortices is however modified by the eflfect of an infinite 
array of image charges across the edges of the finite array of length L and width W . 
Hence, we calculate the net force between the two vortices taking into account the 
effect of the image charges h 

Due to finiteness of the array in the x direction there exists a one-dimensional 
infinite array of images of the charge placed at 0 at points nW n = — oo, • • • 0, • • • oo. 
Since the set of image charges are symmetric about the y axis, the force exerted on the 
test charge placed at P at a distance yij from 0, has only a vertical component given 
by 


*We neglect the effect of the linear defect in the array 


on the vortex-vortex interaction. 
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It is interesting to note that as the above expression is dependent only on y.y as long 
as both the charges are placed along the y axis. 

Having found the effect of image charges due to finiteness of the array in the x 
direction, we proceed to calculate the effect of image charges due to finiteness in the 
y direction. In Fig. 5.4 we show the positions, the descriptive scheme for naming the 
charges and their signs. All distances are measured from the origin. We divide the 
region into two halves separated by the origin and the charge on the lower (upper) 
half is called y^i (2/mu)- Moreover we note that each image charge is in actuality 
a representative of a set of line charge perpendicular to the x = 0 line due to the 
finiteness of the array in the x direction as noted above. The general recursion relation 
that holds between the upper and the lower half sequences of image charges are: 


Vnl = |2/(n-l)I,| 4- L (5.2) 

Iz/nul = y{n-l)l + L 

Vnu — ~|2/(n-l)/l 

yu = L + yi 
yiu=L- yi 


The terms arising out of •4-ve charges in the upper half are from the previous analysis 
of the form cosech{2mL + Ay) while those arising arising out of -ve terms are 
i:~=o cosech([2m+l]i+y), where Ay = y 2 ~yi and Y = y^ +y 2 . A similar analysis for 
the lower half shows that the sum from +ve charges gives cosech([2m + 1]L - Y) 

while those arising from -ve charges give Em=i cosech(2mL - Ay). Incorporation of 
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Position Sequence Chaxge 


-2L A yi y2u 


-L - yi yiu + 


-LI2 

0 ^ 

L/2 


yi you — yoi — 


L-yi yii + 


2L + 2/1 y2i 


X — Q 

Figure 5.4: The figure shows the positions of the image charges in the y-directmn due 
to a -ve charge placed at yi- The positions as measured from the ongiu (^0), the 
sequence for naming the charges and their signs are given on the right. Note that each 
image charge in actuality corresponds to a series of charges in t e x irec ion wi 

same y coordinate. 
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the self-interaction term in the above formalism is straight-forward and is done by 
putting yi = y 2 and hence Y = 2^2, Ay = 0. Taking care that a charge cannot 
have self— interaction with itself we finally arrive at resultant force due to vortex- vortex 
interactions as 


2‘kEj 


TT 


W 


^ cosech ~{2mL + Ay) - cosech + 1]T + 

lm=0 ^ 

■f cosech ^([2m + l]L - F) - (1 - S^o) cosech ^(2^1 - Ay) 
— cosech — ([2Tn -f- IjT -f 2y2) -I- cosech —([2m + 1]L — 2y2) 


(5.3) 


Knowing the force due to such a configuration one can deduce the “current-equivalent” 
by using the relationship F = (/(?!>o)/a where a is the dimension of the plaquette. For 
the general case of two vortices present in the CC at distances y( and y^ one can write 
this force in a more compact notation as 


FiAvly']) 


^ cosech{y'ij) + ^ j {cosech(2A:7' -f jylj) 

-t- cosech[(2fc — l)r — jY-j] + cosech[(2it — l)r — 2yy,-)] } 


(5.4) 


where r = ttLIW, y[ = TcyjW, Y^j = y,- + y'-, and Fij is the force in units of 27rEj/a. 
For either vortex, it now depends upon the direction of movement and is larger for 
moving the vortices closer. 

For multiple vortices present in the CC, the net force Fi on the vortex due to 
all other vortices is given by Fi = Hence, it is given by 


ait = Fi + ip (5.5) 

Here, Fi is taken to be negative (positive) when it moves the vortex in a direction 
.parallel (antiparallel) to that of Fi. A plot of it against Fi is shown in Fig. 5.5 for a 
variable number of pinned vortices and array sizes. For the two— vortex case the force 
Fi is varied by changing yij. This can be done by applying the local current drive id 
to either vortex until the desired yij is achieved. For the multiple vortex case, a range 
of Fi values is obtained by simply measuring it for each vortex in the CC. The fit to 
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Figure 5.5: The threshold current required to depin a vortex plotted as a function of 
the force due to VVI on the vortex. The solid line is a linear least square fit to all the 
measured points using Eq. (5.5) with ip = 0.106. 


Eq. (5.5) is excellent and yields an ip = 0.106, in agreement with the value for a single 
vortex. We conclude that pairwise interactions taking into account images of vortices 
across the edges of the array are adequate to describe all the WI relevant here. 


5.4 Equation Of Motion Of Vortex 

We can now analyze the motion of the vortex— street in the CC. The vortex farthest 
from the defect in array B invariably lies in a region free from fields due to the 
applied in the accelerator region. The equation of motion for a vortex moving under 
an ill ^ square array [137], modified to include the WI, can be written as 
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wii6re nx — Tr/?^ and t] — tt dirnensionlcss vortex rnass and viscosity, respectively Tlie 
coordinate of the moving vortex is y and time is scaled in ^ . The current equivalent, 

of the force on the vortex due to the other vortices, located at positions yj in the 
CC, is given by = Ej F{y, Vi)- should be noted that a single massless vortex, i.e. 
without the VVI and m = 0, in a zero external field = 0), Eq.(5.6) yields y = 0, 
implying no motion. 

In Fig. 5.6 we show the measured position of the vortex in the CC as a function of 
time in array 5 for T = 0 and = 0, 10 and 100. The vortex is pushed by the 12 to 
13 vortices behind it (one gets created at the tip of the defect before the vortex under 
study leaves the array). Also shown are the results of integrating Eq. (5.6) with and 
without the mass term, 7r/?c, for i^xt = 0 and ip = 0.106. We note that m = 0 gives an 
excellent fit to the /3c = 0 and 10 data. For /3c = 100, the m = 0 curve is much closer 
to the experimental points than the one for m = 7r/3c. The data and the m = 0 fit for 
/3c = 100 can be made to agree with an enhancement in viscosity by a factor ~ 3[152, 
153]. We conclude that (i) the VVI are sufficient to explain the motion of vortices in 
the field-free detector region and (ii) the effective vortex-mass is insignificant for the 
studied range of /3c. Our results thus confirm the findings of Hagenaars et ai!.[137] that 
vortices have a negligible mass for /3c < 35. 

A few additional observations are in order. Firstly, while a vortex is moving through 
the field-free region, i^,„ stays just above the threshold ip, as can be deduced from the 
step-like structure in Fig. 5.6. Secondly, the separation between vortices in the CC 
is quite non-uniform, e.g., it varies from 4 to 19 lattice spacings in array B as one 
moves from the defect tip to the edge of the array. It is worth pointing out that in a 
triangular lattice, where the pinning potential is ~ 5 times smaller than in a square 
array, vortices would flow even if they were approximately 50 lattice spacings apart 

^Hagenaaxs et a/.[137] have found that rj is quite nonlinear and can be fitted well by the function 
T] = Al[l + By] where A and S are constants weakly dependent upon /3c. This fit is, however, over a 
very large range of drive current whereas our drive current equivalent of the vortex-vortex interaction 
is very close to ip. The parameters A and B are dependent upon the range of drive current used in 
fits. It IS not clear how well the above form of viscosity fits to the very narrow range of drive current 
applicable to our studies. Hence, we use y = ir. 



nnA PTER 5. BALLISTIC VOBTTCEf^ 


99 



Figure 5.6: A plot of the position of vortex versus time. The origin of time is axbitra^. 
The frequency u, = c. for = 0 and w = u,, for /3 > 0. ^ Imes^are to 

obtained from integrating the equation of motion pven by Eq. (5.6) 

L = 0.106 with and without the mass term. The ^ values are marked next to each of 
the fitted curves and the suffix m denotes inclusion of the mass term. 
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with Fij oc llvij- This is consistent with the estimate of ZFOM that there is at most 
one vortex at a given time in the channel region which is 20 lattice spacings long. In 
any case, as we have shown the vortices in the detector region pile up until the vortices 
are close enough to overcome the pinning potential. Since the VVI is a mutual force, 
there may arise a confusion as to how such a term can lead to a overall drift velocity 
in a particular direction. This seemingly paradoxical viewpoint can be brought in 
agreement by noting that the defect (which acts like a gun shooting out vortices) is 
constantly pushing in new vortices causing a net force in a particular direction. We 
have noted that if fhe current is stopped the vortex street expands due to the mutual 
repulsion between vortices. 

For the sake of completeness, we have also simulated H-shaped arrays with vortices 
produced by an external magnetic field and used bus bars as in the experiments of 
ZFOM. We find the same results qualitatively. Good quality numerical data would, 
however, require prohibitive amounts of computer time This is unwarranted as 
all the essential features of the experiments can be interpreted in the context of the 
experimentally realizable geometry we have used. 


5.5 Discussion and Conclusion 

In a capacitive triangular array, the effect of oscillations of a vortex inside a single 
plaquette is seen to be small ; assuming a typical distance of 20 lattice spacings between 
two vortices, the maximum amount of tilt from the line of alignment is tan-^(l/20). 
The busbars also help in aligning the vortices which as a result do not “fan” out. 
Furthermore, since there is an order of magnitude difference between temperature 
at which vortices diffuse in square arrays as compared to that in triangular lattices, 

experiments probing such effects in square arrays, must be carried out at a much lower 
temperature. 


®The present sunulations were done on a Hewlett-Packard series 9000 model 735 workstation. 
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To conclude, we have shown that VVI can overcome the vortex-pinning potential in 
J JAs and cause the vortices to move through a field-free region. The experimental ob- 
servations of ZFOM can thus be understood in terms of WI alone without recourse to 
ballistic motion of vortices. Our findings are consistent with the numerical simulations 
which have failed to observe ballistic motion for a single vortex. More experiments, 
preferably with a single vortex, are needed to confirm if ballistic motion of vortices 
indeed exists for some other range of parameters. Further theoretical insights as well 
as new experiments are also required to explain the absence of a vortex— mass for the 
range of /3c studied and better understand the dynamics of vortices in JJA. 


f 



Chapter 6 

Results and Discussions 


We now summarize and review our results and discuss their ramifications. 

6.1 Results 

In listing these we proceed chapter-wise. 

Chap. 2 deals with fast simulational algorithms for square JJAs described using 
the RSJ model. We have shown that the complexity of the algorithms can be reduced 
from 0{N‘^) to ^(iVln N) at each integration time step for a number of experimentally 
studied cases. We have done this by noting that there exists an analogy between our 
differential equations and those of electrostatics, that a knowledge of the corresponding 
eigenvectors and eigenvalues of the underl 5 dng differential operators when subject to the 
correct boundary conditions determines the Green’s function (G) and that gauge fixing 
removes the zero mode invariably present in G. For finite arrays we have determined 
the solutions by using a combination of the degenerate eigen vectors which possesses 
the symmetry of the underlying lattice. In all cases we have specifically identified the 
fast transform to be used to integrate the equations of motion. 

• The Periodic Square Lattice: In this case the eigen vectors have been shown to be 
~ exp(ik ‘ fj with A:,- — 2n,-7r/L, i = x, y n,- = 0, 1 • • • L,- — 1 while the eigenvalues 
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are 4 2 cos 2 cos ky. A 2D Fourier transform has been shown to provide a 

fast algorithm. 

• A Finite Square Lattice : The extra Neumann— like boundary conditions arising 
from the finiteness of the array lead us to a solution of the form ~ cos kxX cos kyp, 
kt = n,-7r/I,, * = s, y n, = 0, 1 • • • L,- — 1 for (a:, y) located at the dual lattice sites. 
The eigenvalues remain the same as above. For and Ly of the form 2^, a 2D 
Cosine Transform, provides a faster method. 

• A Finite Square Lattice with Single Bus Bar; By observing that the vanishing 

of wavefunction at the bus bar automatically removes the zero eigenvalue, we 
have been shown that 4^(x, y) ~ sinkxXcoskyy with k^ = {nxT^)/{Lx + l/2) = 

0,1 ••• Lx — 1 and ky as above. These eigenfunctions, however, make the use of 
fast algorithms of the form described above impossible. Consequently, we have 
introduced a variant of the above method to make the process faster. This uses 
cosine transforms in one direction {y) and cyclically reduces the non-diagonal G 
matrix in the k space to a tridiagonal form before solving it in the other. The 
essence of the procedure is described in Fig. 2.5. 

• A Finite Square Lattice with Double Bus Bars: In this case too the busbars in the 
x-direction modify the eigen functions to ~ sin kxX cos kyy with ki = n,-7r/L,- x = 
0, 1, • • • La; — 1, y = 1/2, 3/2, . . .Ly — 1/2 while the eigenvalues remain unchanged. 
It then follows that a fast algorithm is possible if Lx + 1 = 2^ and Lx — 2“^. 

In the case of missing bonds, the above procedure is rendered inapplicable, since 
these introduce internal boundaxies. An alternative route, used by us, relies on isolating 
the effects of bond removal into a separate matrix h and proceeding perturbatively. By 
adapting a procedure commonly used for the case of disordered systems, we have shown 
that the resulting series (see Eq.( 2.49)) can be summed exactly leads to a algorithm 
with complexity 0{NhiN + N). The extension of the above procedure to the case of 
n missing bonds in the form of a linear defect, has been shown to yield a process with 
complexity 0(iV"ln + nN). 
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Chap. 3 extends the above analysis to the case of triangular arrays. We have used 
the Dirac bra— ket notation to describe the JJA and in the process made the mathemat- 
ical analysis much more elegant and transparent. It has been observed that the very 
nature of the triangular lattice makes the corresponding eigenfunctions non— separable 
in the spatial directions. We have noted that the eigen values of the matrix Gg ^ possess 
the 12— fold symmetry of the hexagonal group. For periodic triangular lattices, we have 
had to additionally incorporate the concept of different kinds of periodicity to specify 
the connectivity of the lattice points at the boundaries. We find that there exist two 
forms of periodicity viz., rectangular and hexagonal. We have evolved fast algorithms 
for each, in addition to one for a finite lattice as discussed below: 

• Rectangularly Periodic Triangular Lattice: The results for this case have yielded 
a solution is similar to that for a periodic square lattice with the exception of the 
eigen values which are now 6 — cosfcj — cosfc2 — cos(fci — ^ 2 ); K = 2n,'7r/L,i = 
0 ,1; Tii = 0,l,---i,- - 1. 

• Hexagonally Periodic Triangular Lattice: Viewing the josephson lattice as a sub- 
set of a 2D-plane on which data is finitely sampled, we have shown that for such 
a lattice the eigenfunctions are ~ exp— i ^(2x1 — X2)(2ni — 712) -f £2:2^2 and 
possess the 12-fold symmetric eigen value 6 — 2cos(2ni — n2)|^ — 2cos(2n2 — 

“ 2cos(ni +712)11 with 0 < tii < 31, 0 < 712 < I. The properties 
(see Sec. 3.3.2) of the wavefunction have been shown to yield a fast algorithm if 
L = 2P. 

• Finite Triangular Lattice: In this case we have shown that that the wavefunctions 
do not completely diagonalise the matrix . However, after certain manipula- 
tions, Gq ^ is seen to reduce to a tridiagonal form which can then be solved with 
increased efficiency. 

Chap. 4 deals with our work on bond diluted and disordered defects. The results 
are as below: 
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• For the case of bond disorder, with grading g < gm = 0-2, we find vortices pinned 
inside the linear defect. The pinned vortices are seen to increase the stability of 
the system against breakdown. 

• Our findings reveal that the barrier height for vortex motion, and hence the 
depinning current, is highly anisotropic and can be greater than for rectangular 
arrays. As a result, we observe that vortices are pinned by a long narrow lattice 
with bond-diluted linear defects of size much less than predicted earlier. These 
pinned vortices are seen to raise the ic of the system (Fig. 4.5). We also show that 
there exists a scaling behaviour of iAicjlT Nxl^ for fixed A'y/n (Fig. 4.4). 
Moreover, for small Nx^ vortices are seen to flow in the direction of the current 
due to an increased gradient of current flow near the tip of the defect, 

• The energy of the system in the steady-state regime has been conveniently di- 
vided into various sectors with each sector being characterised by the number of 
pinned vortices present in the array. The energy decreases continuously within 
each sector and registers a discontinuous downward kink while crossing from one 
sector to a higher one. Hysteresis in the E vs iext plane is observed. This occurs 
because the pinned vortices do not escape from the array, even when the external 
current is reduced to zero. 

• The periodic regime of the array reveals the presence of waves that carry the 
vortices. The wavelength of the wave is seen to increase as one moves towards 
the edge of the array. A snapshot of the phases also reveals a “wake” that moves 
ahead of the vortex and provides energy to propel it forward. 

• An analysis of other defects patterns affirms the analogy between current flows 
in JJAs and hydrodynamic flows. Since the curvature of flow lines decides the 
point of breakdown, effective streamlining of a defect is seen to increase v- 

We have investigated in Chap. 5 the flow of vortices predicted to occur ballistically 
in underdamped JJAs at finite temperature. Our findings axe summarised below. 



106 


• We have shown that a linear defect provides a suitable geometry to experimen- 
tally and numerically study ballistic vortex motion since it produces a naturally 
collimated beam of vortices. 


• It is observed that although a single vortex does not move ballistically, a vortex 
street can move through a field-depleted zone even in an overdamped array. 

• At finite temperatures we find that the vortices diffuse so that only a fraction 
of those produced by the defect reach the detector. The discrepancy, between 
our simulations and experiments in the literature, in the temperature at which 
such diffusion becomes prominent, is argued to be due to a difference between 
the pinning potentials of triangular and square arrays. 

• We have found that the current required to depin a vortex in the presence of 
other vortices is ~ O.lio as found for the infinite array limit. To arrive at this 
value, we take into account the complete vortex-vortex interactions inclusive of 
the effects of the array of image charges present across the edges of the lattice 
(Eq.( 5.4)). 

• We determine the displacement of a free vortex as a function of time by including 
the current equivalent of the WI in the vortex equation of motion (Eq.( 5.6)). 
The agreement between the numerical and observed data is startling (Fig. 5.6). 
Our curve further reveals that for /?c ^ 30 the mass of the vortex is zero. 

6.2 Discussion and Scope For Further Study 

We have shown that algorithms for simulating a number of experimental JJAs e.g. 
those with bus— bars and/ or defects as also those subject to magnetic fields, noise, 
and/or current drives, all involve Green’s functions specific to the underlying lattice 
in question. Furthermore, the Green’s function of interest can straight— forwardly be 
constructed once the eigen-basis of the corresponding connectivity matrix or discrete 
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Laplacian has been determined. We have provided a much more transparent rederiva- 
tion of the Eikmans-Himbergen algorithm using this method. The extensions of this 
algorithms to the case of 3D arrays is easy since the corresponding eigen-vectors are 
specified to be either exp(ik ■ f) or cos[kxx) depending on whether the boundaries 
are periodic or free boundaries respectively. The eigen-values are now changed to 
6 — 2coskx — 2 cos ky — 2 cosk^. In this case, one has to obviously use 3-D transforms. 

We have also provided a general method for adding and removing bonds from a 
configuration whose Green’s function is known. This can be used to create a variety 
of defects whose properties can then be studied. In this thesis we have only explored 
the possibility of bond-removal. It would be interesting to determine the effect of 
bond-additions especially since these would change the pinning potential for vortices. 
Random defect configurations which can also be studied by this method, represent yet 
other systems of interest. 

It is also worth noting that many of the expressions we have derived, e.g. those 
for the Green’s functions, can be written in terms of continuous functions, which are, 
of course, to be evaluated or sampled only at points belonging to the array. Now if 
we turn this observation around and think of the discrete lattice as being produced by 
the sampling of a continuous 2-D waveform, we make contact with a long-standing 
problem in electrical engineering, viz. that of the digitisation and subsequent recovery 
of band-width limited analog wave-forms. We can thus make a useful correspondence 
between the algorithms presented here and those used in digital signal processing of 

data. 

We have used this analogy, in part, to extend the scope of our algorithms to tri- 
angular arrays. A case in point is that of the hexagonally periodic/finite triangu- 
lar lattice which in either cases admits eigenvectors akin to those derived arrived by 
Merserau[133]. The finite case is rather more difficult for the triangular lattice than it 
is for the square one. This is because the creation of a boundary in the former, requires 
the removal of two bonds placed obliquely with respect to each other, rather than one. 
This case therefore requires a new approach and this has been developed in this thesis. 
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In particular we have broken up the inverse laplacian Gq into a periodic (P) and a 
boundary (S) matrix and have shown that for a specific choice of the basis both B 
and P — B — take on forms which are amenable to diagonalisation in linear order. 
This has allowed us to present a computationally fast algorithm for the finite case. 

It is also been pointed out that for a hexagonally periodic lattice, equivalent points 
in different tiles (see Ref. [133]) can be exchanged to create a new fundamental pe- 
riod. As a result these periods come in several shapes of which Fig. (3.1) is but one 
example. The corresponding eigenvalues and eigenfunctions of Gq^, however, remain 
unaltered and the O(NlnN) algorithm continues to be equally applicable to each of 
them. These shapes can serve as starting points for creating several different finite 
lattices. The various types of boundaries which result are especially significant in view 
of the interesting question: For which class of boundaries does the matrix B commute 
with P? Indeed if there are cases for which [P, B] = 0, they will automatically lead to 
0(N In N) algorithms. 

The case of the triangular lattice opens up a set of questions which can form the 
foci of future investigations. For example, it would be interesting to determine the 
differences in the ground state of a hexagonally periodic and a rectangularly periodic 
triangular array in the presence of a magnetic field. These differences, if any, would 
reflect itself in a change of Tbkt- Furthermore, it would be useful to find if accidental 
degeneracy discovered in the case of rectangularly periodic arrays[139] with / = 1/3 
and / = 1/4 occurs in the hexagonally periodic lattice too. Further investigations on 
the presence of various classes of Shapiro steps[14l], the role of positionally disordered 
defects in dictating the dynamics of arrays[74] and more importantly the detection of 
ballistic motion in arrays of larger sizes can be attempted. 

For the case of defects, several directions for further research are possible. Firstly, 
an understanding of the steady— state regime from a fixed point analysis of the evo- 
lution equations of the system is desirable. This can be achieved relatively easily for 
small arrays but becomes formidable as the system size increases. For an A'* = 4 array, 
a knowledge of the values of the external current and the phase at (—1/2, Ny + 1/2) 
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(Fig. 4.1) provides a complete solution via the application of current conservation at 
each node. Hence is a function of only one variable. A plot of lext versus the 
concerned variable provides solutions which consist of purely superconducting flow. A 
careful study, further, reveals the that unstable solutions are also present. The break- 
down of steady-state flow (4) is characterised by a marginally stable solution. The 
challenge is to generalize this procedure to larger arrays. In the case of pinned vortices, 
the position of the vortex constitutes a further input into a complete solution. Secondly, 
it is of interest to characterise the spin waves seen in the periodic regime, in terms of 
their dispersion relation. The subtlety lies in the fact that there exists a “frozen-in 
wave in the steady because of which k 0 as cv —>■ 0. A similar deviation from con- 
ventional spin wave behaviour exists in the presence of pinned vortices. Thirdly, an 
in depth study of the symmetry-breaking of the solutions d{x,y) = —9{—x,y) in the 
chaotic regime is in order. In fact, we feel that the breaking of this symmetry, is a 
much better and faster signal for detecting chaos in JJAs than the conventional power 
spectrum method. Phenomena like bunching and debunching of several vortices in the 

column-switched regime also merit deeper study. 

The idea of a “vortex-street” moving ‘ballistically’ through a field depleted zone has 
been introduced in Chap. 5. It has been shown by us and other authors independently 
that the mass of a vortex is negligibly small in contradiction to theoretical calculations. 
It would be interesting to see whether the interactions with the image charges due to 
the finiteness of the array, actually renormalise the inertial mass of vortices to zero. 
Since the exact' position of the vortex inside a plaquette is unknown, one has to follow 
a boolean kind of logic whereby each vortex is either midway between the plaquette 
or absent. This is important if one tries to find the vortex motion of the second last 
vortex (instead of the last). Here any error in the position of the vortex is magnified 
by the VVI and makes the comparison between theoretical and numerical results much 
more difficult. It would hence be of interest if one can determine the position of a 
vortex inside a plaquette accurately. 
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